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Those who aspire to inspire, perspire 
Ever foiled, embroiled in unending toil 
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In gleaning the meaning, whilst ever daydreaming 
Of the problem they're solving, dissolving 
And a time when the rhyme rings sublime 
Lest their sighs, hide the prize, from their eyes. 



Preface 


My undergraduate maths course on Number Theory was given by a 
particularly enthusiastic lecturer. There were about a hundred students on 
the course and the lecture theatre had three large blackboards across the 
front. This became the setting for a most memorable maths lecture. Starting 
on the extreme left of the first board, the lecturer invited the class to call out 
single numbers. As we did so, he wrote them onto the boards as one long 
integer. When we tired of calling out digits, the lecturer stood back and 
surveyed the resulting integer. 

Then the magic began. Without stopping for breath, and without any written 
working, the lecturer wrote down the prime decomposition for the enormous 
number. By the time the problem was reduced to the final pair of primes, 
there was only one student, somewhere at the back of the class, who had 
managed to keep pace. 

When tutoring maths I noticed that factorising seems to be a process which 
variously frightens, fascinates, or delights my students. The simple 
systematic rules which I have long taken for granted are a novelty to them. It 
is rather less than obvious to a beginner that the square root of the number 
sets an upper bound in the search for factors. When I tell them that they 
need only restrict their attention to prime divisors is as though they had been 
given a box of chocolates. 

Inspired by the responses of my students, my own curiosity became the 
motivation for this work. The following account is offered as my exploration 
of the subject of factorisation insofar as it applies to integers. I have 
described the application of simple algorithms which seek to accelerate the 
solution of the problem of factorisation, and have described some physical 
analogues. 

The work is not intended to be exhaustive. For those interested in this 
particular mathematical challenge, there are many avenues for development 
of the procedures described. The mathematics employed should be 
understandable by able senior high school pupils and amateur 
mathematicians. It should certainly be well within the grasp of first year 
university students. The contents may be of interest to those studying 
Computing Science, Mathematics, Engineering or Physics. 

J W Cahill 
September 2007 
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Introduction 


Initially assume that a simplistic approach is used in the search for a 
factorisation of a large number. Then starting with the number 2, each 
successive prime number in turn will be divided into the target number in 
order to establish whether or not it divides exactly. 

There are, of course, some well known shortcuts. If the target number ends 
in 2, 4, 6, or 8, one of its prime factors is 2. If it ends in 5, one of its prime 
factors is 5. In either case, the problem is trivially reduced. If the target 
number ends in 0, it has prime factors of 2 and 5. A further trivial test which 
can be carried out is to add all the digits of the target number together. If the 
sum is divisible by 3, then the target number itself is divisible by 3 and again 
the problem is simplified. 

The same process is applied to the quotient obtained from the first pass and 
the exercise repeated. Provided the target number was highly composite, it is 
possible to make rapid inroads to the factorisation. 

If the target numbers are large and less composite, the issue of factorisation 
takes on a rather different character. The question of just how long it might 
take to factorise a large number is not long in raising its head, and recourse 
is made to computers. 

In the worst case, the number turns out to be prime. Should that happen, 
then attempted divisions range up to the nearest prime less than the square 
root of the target number. The difficulty of factorising large numbers is well 
known. Not only does the number of potential divisors increase, the 
complexity of the division itself grows in tandem. Consequently, the time 
taken by a computer to carry out a division increases significantly as the size 
of the numbers increases. 

In this study, a sieve of Eratosthenes was incorporated in one of the 
programs which was used for generation of primes. The sieve removed all 
potential numbers divisible by primes up to 29, but the effectiveness of such 
a sieve rapidly diminishes. There does not appear to be any obvious method 
to achieve a dramatic reduction in the number of possible divisors. The 
generation of primes in the range 1 to 10 7 identified approximately 1.6 x 10 6 
possibles, roughly 16% of all the numbers in the range. When it is 
considered that ignoring even numbers eliminates 50%, of all possibles the 
inclusion of an additional 9 prime numbers only eliminated about 70% of 
the remaining possibles. Of the potential primes identified, about 6.6 x 10 5 
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were found to be prime. So in the range used, roughly half of the integers 
identified as possible primes turned out to be prime. 

Division of a number to be factorised can be carried out regardless of 
whether the divisor is prime or not. If the divisor is large, this approach may 
be less arduous than ensuring that the divisor is a prime before attempting a 
division of the target number. 
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A First Approach 


Prime numbers are not uniformly distributed, there being a greater density of 
primes at the low numbers. This favours algorithms organised to range 
upwards from unity. For a target number, T, selected at random, the average 
factorisation search would tend to require less than Jr 12 attempts using 
such algorithms. A correspondingly greater number of operations would be 
required for algorithms which range downwards from V T. The basic 
algorithms as described here generally range downwards from Jr. 
However, the algorithms may be started from any convenient point or could 
be modified to range upwards. Consequently no distinction is made between 
the two approaches and in terms of average performance, a factorisation 
could be expected by the time the divisor had reached JT12. If it is assumed 
that the divisors are generated by successive increments and that about 75% 
of the divisors so created are discarded, the average number of operations 
required for a search carried out by direct division would be: 

G D /2 ~ ~^2T (Increments) + ^ ( Divisions) 


In the general case of large integers, multiplications and divisions carried 
out on conventional computers require a number of clock cycles. The 
number of operations taken to complete a multiplication or division 
increases as the size of T increases. It would be attractive to eliminate these 
processes as far as possible from the factorisation operation. 


When the target number is divided by some divisor D, the division 
algorithm produces a quotient Q and a remainder R, thus: 

T=DQ+R 


In attempting to avoid the division process, the question arises : 

C \ 

What is the relationship between the 
quotients and remainders from adjacent 
integer divisors? 

v_ ) 


The algorithms described seek to exploit the relationships between divisors, 
quotients and remainders. 
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A Basic Decrementing Algorithm 

Let T be the target number to be factorised. 

Determine the value of the square root of the target number as a starting 
point. Round down to the nearest integer. This is taken as the divisor “D”. 

The target number is divided by D. This gives a quotient “Q” and a 
remainder “R”. This is the only full length division operation required. 

D is loaded into a register which is to be decremented. With each 
decrement, the value in the quotient register is added to the remainder 
register (Q + R —> R). The result is compared to the contents D of the 
divisor register. If it is greater, then the divisor register value is subtracted 
from the value in the remainder register to give a new remainder (R - D —> 
R) and the quotient register is incremented by one (Q + 1 —> Q). The 
operation is repeated, incrementing the quotient register each time until the 
remainder register contains a value which is less than the value in the divisor 
register. 

If, at any stage, the remainder register = 0, then the divisor register and the 
quotient register contain a factorisation of the target number. This algorithm 
will ultimately produce factors. The process is illustrated in the flowchart of 
Figure 1. 

The starting point is entirely arbitrary. The square root is a natural starting 
point since it is a geometrical “Mid point” for the search, and at least one 
factor must lie between zero and the square root. However, all that is 
required is a selection of some value D from which Q and R are determined. 
Once these values are loaded, the computation can proceed uninterrupted 
until it is terminated when R = 0. Termination is guaranteed to occur for 
some numbers D and Q. 

How the algorithm works 

Taking the starting point to be a divisor D which is the nearest integer less 
than the square root of the target number, the quotient Q will be at least 
equal to the divisor at the start of the algorithm. A change of unity in the 
divisor register corresponds to a total change, equal to the quotient register 
contents, in the product of the divisor and quotient. Consequently the 
remainder must be changed by this amount. The remainder is then greater 
than the divisor register contents. This is because the quotient is at least as 
large as the divisor and the divisor could therefore have been divided into 
the target number at least one more time. Hence the quotient register has to 
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be incremented by at least unity and the remainder register reduced by a 
corresponding multiple of the value of the divisor. Since the divisor register 
is being decremented, and the quotient register incremented, the value in the 
quotient register is always the greater. 



The basic decrementing algorithm 
Figure 1. 


Assume that both registers initially contained values close to the square root 
of the target number. As the search continues, there will come a time when 
the divisor register contents will need to be subtracted from the remainder 
register numerous times before the remainder value will become less than 
the divisor register value. The contents of the quotient register will grow in 
proportion. 

Speed of the algorithm 

Assume that the factorisation starts at the square root of the target number. 
In the worst case, the target number is prime and the quotient turns out to be 
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the target number itself. Thus taking the target number as T, the number of 
loop repetitions in the worst case is: 

k = T- JT 

The divisor is decremented from the square root of the target number down 
to unity. The remainder undergoes the same number of additions as the 
number of decrements in the divisor. It also undergoes k subtractions while 
the quotient undergoes k additions. Thus in the worst case, the total number 
of add and subtract operations is: 

G wom * 2(T-jT) + 2jT 
G worst ~ 2 T 

The speed of the algorithm is not constant. The performance is best when 
the divisor and quotient are of similar size. As the divisor becomes very 
small, the algorithm slows down dramatically, with the worst case being 
when the divisor is set equal to unity. This final step will require 
approximately T/2 loops with corresponding additions and subtractions 
amounting to T operations. However, when the algorithm has been 
completed for the divisor = 2 and has failed, the number is known to be 
prime. Consequently it is not necessary to carry out the final stage of the 
search algorithm. Similarly a test on the final bit of the target number 
determines whether or not it is even, eliminating the need for the 
penultimate step as well. Ignoring the single decrement in the divisor and 
the corresponding single addition to the remainder, the number of additions 
and subtractions required between a divisor of 3 and completion of the test 
with the divisor set to 2 is: 



The test for a divisor = 3 can be obviated by adding the digits of the target 
number and dividing the result by 3. Likewise the final digit of the target 
number can be tested for divisibility by 5. This leaves 7 as the lowest prime 
which needs to be tested explicitly by the algorithm. For small values of n, 
the total saving in computation is approximately equal to the sum of these 
economies: 
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" n..(n -\) : 


• 1 ~ 27 ’{( 1 2 ) + C 2 3 ) + •••• + [ («-2) (n-1) 


So for the case where the lowest number tested is 7, the saving is: 


Consequently, the number of arithmetic operations required for the 
algorithm to test all divisors of practical interest is: 


In practice, the tests for 2, 3 and 5 are just as well carried out by direct 
division as by bit testing since they are one off tests. Furthermore, the 
computational savings described above apply under conditions in which 
direct division would be preferable in any event. Consequently, in arriving 
at an estimate of the number of operations required for a typical 
factorisation by the algorithm using sequential computer operations, these 
savings are ignored. As before, the average range which must be covered 

IT 

before a factorisation is found is taken to be —. It is readily shown that the 
number of operations required is then given by: 

[T 

Gdi 2~2 \- ’+2 JT = 3 JT 

The algorithm as described decrements the divisor register by unity each 
time. The only divisors of interest are odd numbers, so having chosen a 
starting point which is an odd number, it is a straightforward matter to 
decrement the divisor in steps of 2. 

The relationship between successive decrements is then given by: 

D j =D j~ 1“ 2 
Qj = Qj-\ +2 


Rj -Rj ~i +2(g/-i ~Dj_ i)+4 
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With this change implemented, the average number of operations to achieve 
a factorisation is reduced to: 


g DI2 


3 JT 
2 


A Preliminary Improvement 

Since the basic decrementing algorithm described slows down markedly at 
small divisors, a simple refinement of the strategy uses direct division of the 
target number for several of the smaller primes and employs the result to 
limit the progress of the basic algorithm. Although in principle it is 
necessary to divide the target number by each of these prime numbers in 
turn, it is also possible to multiply sets of prime numbers together and divide 
the target number by the resulting composite number. The remainder from 
this division is a number which is smaller than the product of the set of 
primes. 

The remainder can in turn be divided by each of the primes from the set. If 
any one of these primes divides the remainder, then it also divides the target 
number. Division of the remainder is a much smaller problem than division 
of the target number itself and hence can be carried out more quickly. 
Implementation of this refinement excluding the primes up to 19, showed a 
noticeable improvement in performance when attempting a factorisation of 
primes which were greater than about 10 10 in size. 

Relative Growth of the Quotient 

If a multiplication was carried out, the product of the divisor and quotient 
would always be close to the target number. Thus: 

T ~ DQ 

If D and Q are approximately equal initially, then following a decrement in 
the algorithm, D will only have to be subtracted once from the remainder to 
ensure that the remainder is again less than D. Likewise Q will only require 
to be incremented by one. 

However, as the remainder grows, it will occasionally be necessary to make 
a further subtraction of D to return the remainder to a value less than D. At 
these points there will therefore be an additional corresponding increment in 

Q. 

As the search progresses, there will come a time when Q is twice the value 
of D. From this point until Q = 3D, the divisor will have to be subtracted 
from the remainder twice, with occasional extra subtractions as above. In 



general, the quotient and divisor are related by: 

Q=fiD+8 

For some integers /i: l <ju< JT and 8 : 0 <S < D. 

Each increment in the multiplier p increases the number of times that the 
divisor must be subtracted from the remainder. The points at which at which 
these loop expansions occur correspond to: 

D = -yt- JT and Q= Jju JT 

Since: 

T*DQ=-^jTj]u JT 

For each value of p it is therefore possible to calculate the number of 
occasions on which the remainder will be equal to or exceed D. These are 
the only points at which a factorisation might exist, with factorisation should 
the remainder = D. For example, between p = 1 and p = 2, the divisor 
decreases by: 

The quotient increases by: 

AQ = (J2 -1)JT 

Assuming the factorisation is not found in this range, the difference in 
magnitude between these two numbers is the number of occasions on which 
the remainder was greater than D. 

AQ-AD= J2 -2+ ~y~ JT JT 

V ^ 
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Example 1. 

Taking T with some value in the region of 10 6 and using a starting point at 
the square root of the target number, we have: 

T ~ 10 6 

Jf ~ 10 3 

A D * 293 and A Q ~ 414 
AQ-AD « 121 


The quotient is incremented once for each decrement of the divisor down to 
about D = 707. There are some 121 extra increments of the quotient 
corresponding to each time the remainder exceeded the divisor. Only at 
these 121 convergences was there a possibility of factorisation. 

Note that since T is odd, the quotient and the divisor, must both be odd. If 
the algorithm was started with an odd divisor and odd quotient, then the first 
convergence cannot be a factorisation. This is because there is an extra 
increment of the quotient at the point of first convergence. Consequently, 
the quotient would become even. However, the target number is itself odd 
and cannot therefore contain any even factor. 

From the first extra increment until the second one, the quotient will remain 
even. There is no possibility of a factorisation until the next extra increment. 
It is evident that for any searches which commence with both quotient and 
divisor odd, every odd extra increment, cannot provide a factorisation. In 
the above example, from 121 extra increments, there are therefore only 
approximately 60 instances in which a factorisation might have occurred. 

From the decrement of the divisor following D = 707, the quotient will need 
to be incremented twice for each decrement in the divisor. This situation 
will exist between p = 2 and p = 3. There will also be additional increments 
of the quotient every time the remainder exceeds the divisor. The analysis 
can be applied for any given value of the multiplier, p thus: 


A D = 


1 

Jf 


1 


Jf 
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aq = (J07TJ) 


AQ-AD 





1 

yc^ry 


JT 


Example 2. 

Using the basic algorithm, the number 1987654323 was reduced to factors 
9777 (=3 x 3259) and 203299 (= 263 x 773). The starting value was the 
square root of the target number rounded down to the nearest integer. This 
was 44583. The initial remainder was 10434. So the reduction took 44583 - 
9777 = 34806 loops to complete. The initial value of u was unity since both 
factors were initially the same. The final value of li was the nearest integer 
less than 203299/9777 = 20. Towards the end of the operation of the 
algorithm, there were therefore 20 repeat subtractions and the same number 
of increments to the quotient at each decrement of the divisor. 

The u5 Improvement to the Algorithm 

Evidently, keeping track of the growth in p has the potential to accelerate 
execution of the algorithm since several repeat subtractions and increments 
could be replaced by a single operation. Furthermore, the remainder need 
not have the quotient added to it, only to be followed by repeated 
subtractions of the divisor. In place of these operations, the difference 
S = ( Q-juD ) may be added to the remainder at each decrement of the divisor 
D. The difference 6 may be stored in one register and the multiplier p in 
another. 

Once this modification has been incorporated, it is no longer necessary to 
retain a quotient register, and the operations associated with it are discarded. 
Once the algorithm has terminated with R = 0, the quotient can be 
calculated. Figure 2 illustrates this modification in flowchart form. 

This algorithm produces a noticeable increase in speed compared to the 
basic form. However, there are additional counters and comparisons 
required and the detail is slightly more complex than illustrated. These 
factors detract from overall performance in a software implementation. 
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General form of the single step ad decrementing algorithm 







The Direct Quadratic Approach 


In the enhanced decrementing algorithm already described, the quotient was 
replaced by a multiple p and a difference 5. Consequently division of the 
target number may be rewritten: 

D(p.D+d) + R = T 
/j.D 2 +SD + R-T=0 

It is usual for the coefficients of a quadratic equation to be known and for 
the roots to be determined accordingly. In this particular application, the 
only predetermined quantity is the target number T. However, a set of 
constraints are applied which enable a solution to the equation to be found. 

T and R are both integers. By definition, p and 5 are also integers. Thus all 
the coefficients are integers. In those cases where p = 1, it follows that the 
only real integral roots of the quadratic equation are factors of the constant 
coefficient. 


Since the only cases of interest are those in which the remainder R = 0, the 
constant coefficient is simply T. Consequently the roots of the quadratic are 
given by: 


D = 


S±Jd 2 +4fiT 


In which all the variables on the right hand side are integers. The condition 
that the roots, D, be rational, requires that the term (5 2 - 4pT) be a perfect 
square. For the roots to be integers, the numerator must also be exactly 
divisible by 2p. The further requirement that the roots are odd integers 
requires that the numerator be an odd multiple of 2p. 

By inspection it can be seen that for all values of T, the formula returns 
values of D = 1 or -1 when p = T and 8 = 0. Similarly, when p = 1 and 5 = 
T-l, the formula gives values of D = 1 or -T. 

A First Quadratic Search Algorithm 

As before, T is tested to ensure that it is an odd number. Then, p is set equal 
to unity. Since the only discriminant values of interest are perfect squares, a 
variable X is defined by: 

A 2 =S 2 +4pT 

Assuming that the search commences at the square root of the target number 
sets 5 equal to zero. This defines Xo by: 
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^0 = ffiiT 

If T is a perfect square, then X 0 is itself a factor and the process ends 
trivially. Otherwise A,o is rounded up to the next largest integer X. The first 
possible value for 5 is then determined from: 

d 2 = A 2 -4juT 

Since X 2 is a perfect square, 5 2 is the amount which has to be added to 4pT 
to make (3 —4juT) a perfect square as required by the problem constraints. 

Furthermore, the value of 5 2 determined above is itself required to be a 
perfect square in order that the quadratic coefficient 5 will be an integral 
value. This provides the necessary test in the algorithm loop. If 5 2 is not a 
perfect square, then X is incremented and a new value of 5 2 established. The 
loop continues until a value of 5 2 which is also a perfect square is 
encountered. 

In practice, the sequence of incrementing X, squaring it and subtracting 4pT 
is replaced by an appropriate increment in 5 2 since: 

(A + l) 2 = A 2 + 2A+ 1 

.-. (2+ \)~ -A 2 = 2/1 + 1 

(A + 2) 2 = /l 2 + 4/1 + 4 

.'. (A + 2) 2 - (2 + l) 2 = 2i + 3 

And in general: 

(A + n)~ - (A + n — l) 2 = 2A + 2n - 1 

••• dl=8 2 n _ l+ (2A + 2n-l) 

The latter term simply requires to be incremented by 2 in each loop of the 
algorithm. Notice that in this algorithm, X commences with a value 
approximately equal to the square root of the target number and increases 
with the execution of each loop. 

Considering the case in which p = 1 throughout, we have from the 
specification of 5, that: 
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0 < 8<T 


0<<5 2 <T 2 
X 2 = S 2 +4T 

4T<X 2 < T 2 +4T 

The worst case is when T is prime. The number of loop repetitions required 
in that case is: 


Nworst — J A max 

N wors t = Jt 2 + 4T-4T 

■ ■ N worst — T 

Compared to direct division, this is apparently a rather poor algorithm in 
terms of the worst case number of loops. However, even in that eventuality 
the loop itself is very simple, principally consisting of additions and 
comparisons. The most onerous calculation is the determination of whether 
or not 5 2 is a perfect square. 

As before, the average number of loops required to find a factorisation is 
that needed to cover the range JT ( Jt 12). For the case p = 1: 

We have that QxD = T and Xq = 4f.iT 

where Q = fiD + S 


thus when D = — ,Q~1 JT 


giving <5 = {2 - y } JT 
Also X 2 =S 2 +4fiT 




so that = 


X-X 0 = {^-2}jT 


Jt_ 

2 


Lambda need not be incremented in steps of unity, and an improvement to 
the algorithm increments X in steps of 2, reducing the average number of 
iterations to AX = JT 14. 
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Initially 5 2 is very small. The time taken to determine whether or not 5 2 is a 
perfect square will be small. Furthermore, in contrast to the decrementing 
algorithm, 5 2 increases by 0(2),) at each execution of the loop. Thus, in the 
early stages at least, the corresponding values of 5 2 increase rapidly, 
skipping values which could not possibly give rise to a factorisation. 
Nevertheless, this growth benefit diminishes as 5 2 becomes large in relation 
to 0(2/.), and the process is further slowed by the increasing time taken to 
determine whether or not 5 2 is a perfect square. 

The problem of testing 5 2 can be addressed in part by eliminating those 
values which end in 2, 3, 7, or 8. No perfect square can end with one of 
these digits. 

Reverse division, with early termination on failure, as described elsewhere 
in this study might also be applied, since in the case of a square root, both 
divisor and quotient must be equal. A forward square root algorithm can 
also be arranged to terminate on failure, signalled by the presence of a 
fractional part in the square root. 

Performance can also be improved by increasing the value of p at the 
appropriate points in the growth of X. These points correspond to equivalent 
junctures in the progress of the decrementing algorithm and occur whenever 
<5 2 > T/(jU + 1). 

Figure 3 illustrates the behaviour of 5 2 as X is increased through the first few 
increments of p. At each increment in p, the loop requires to be initialised. 
This involves calculating a fresh value for X 0 along with other variables. 
These calculations are required more frequently as the value of p increases 
detracting from overall performance. Conversely, it is evident that the 
progress of the algorithm is much more rapid immediately following each 
increment in p. There is therefore a corresponding reduction in the number 
of loop repetitions. 

The graph of Figure 4 depicts successive differences in S as X is incremented 
in steps of 1. The actual differences depend on the value of the target 
number, and the savings are therefore greater for larger numbers. 
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Graph illustrating typical growth behaviour in 5 squared as X is incremented. The initial rapid growth 
declines as the value of 8increases. The growth rate is restored with each increment in p. 

Figure 3. 



Difference between successive values of 6for the quadratic search algorithm with p= 1 
and unit increments in X. The step size diminishes as 5increases. 

Figure 4. 
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It can be seen from Figure 4 that even with a number of the order of 10 9 , it 
only takes about 100 iterations to reduce the step size in 5 below a value of 
20. Nevertheless the value of § has been increased by 1200 in the course of 
the 100 iterations. This is a useful saving in comparison to either direct 
division or the decrementing algorithm. Although the curve is not 
continuous as a result of the integer nature of the problem, it is worth noting 
its smoothness. 

The test for 5 2 being a perfect square may also be addressed in another 
manner. Having determined some value for 5 2 using X, it is possible to 
identify r|, the nearest integer to, but less than 5. By considerations identical 
to those applied to X we have: 

(tj + s)~ - 0/ + s - 1 ) 2 = 2t] + 2s - 1 

'7.? = '7s-i+( 2, 7 + 2s-1) 

The equation of interest then becomes: 

'7? = <>n 

Where >7? = + (2q + 2 s- 1) 

and 8 2 = (5 2 _ j + (22 + 2 n - 1) 

In its simplest form the algorithm proceeds by incrementing s and n in turn: 

if >7? > Sr, 

Then n <- n + 1 
Else s <- s + 1 
Until r\r = d}, 

The explicit test for S 2 being a perfect square has been replaced by a direct 
comparison to some number q 2 which is known to be a perfect square. The 

value of 5 2 is, however dependent on X and so is not necessarily a perfect 

square. The penalty in terms of computation for this basic comparison form 
is that ,v assumes every integer value in turn. Recall that r) 2 = A 2 -4uT . 
Initially 8, and hence tj <c X. Therefore a large number of increments in s 
will be required for each increment of n. 

Accordingly, this basic form of the comparison algorithm becomes 
appropriate when the difference between successive values of 5 becomes 
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small. Arbitrary step sizes in s and hence r| can be accommodated. For any 
number r e H, we have: 

(>/ + r ) 2 = 17 2 + 2 ip-+ n 

(;; + r ) 2 - 17 2 = 2ijr + >~ 

The value of r 2 remains unchanged for significant numbers of loops, but the 
term 2r|r has to be calculated by multiplication. This may be simplified, 
since for any keU: 

{t /0 + (k+ l)r } 2 -(> 7 o + kr) 2 = iiQ + 2ki]Qr + 2ijQr + k^r- 
+ (2k + l)r 2 - {i/g + 2t)Qkr + k~r ~} 

= 2)]Qr+(2k+ 1 )r 2 

.'. {tj o + (k + l)r} 2 = (>iq + kr) 2 + {2 t/pr + ( 2 ^+ l)r 2 } 

The incrementing value is based on some initial r | 0 and proceeds by adding 
2r 2 to its value at each iteration. Naturally this still requires initialisation 
whenever the step size r between successive values of 5 changes. Referring 
again to Figure 4, it can be seen that the region of the curve in which the 
step size in 5 varies slowly, extends over a significant range. 
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The Behaviour of the Remainder and a Second Quadratic 
Algorithm 

Considering once more the decrementing algorithm, it was shown in 
Example 1, that there were only 121 convergences of remainder and divisor 
with about 60 possible factorisation points out of nearly 300 decrements. As 
each step of the algorithm was carried out, there was a new value of 
remainder. Figure 5 is a graph showing the typical behaviour of the 
remainder as each new divisor in turn is selected. 



200 180 160 140 120 100 80 60 40 20 0 


Divisor 

Behaviour of the remainder as the divisor is decremented in steps of 1. 

Figure 5. 


The variations in the remainder appear to be almost unpredictable. The 
search terminates when the remainder = 0. The divisor is constantly 
changing, and because of this it is convenient to consider the margin 
between the remainder and the divisor. Given the results of a division into 
the target number, the initial difference between the divisor and quotient can 
be determined. 

<^o = Qo~ D o 

The difference between the remainder and the divisor is then the margin: 

mo = (D 0 -R 0 ) 
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During a cycle of the algorithm, the divisor is decremented by unity and the 
remainder increases: 


R - 0 Ro + e o ) 

*!«-(*'-(Z) 0 -1)) 

^1 = R o+(Qo~ D o) + 1 

R l = R 0 +S 0 + l 

However, since the divisor is now one less than it was, the margin has 
become: 

m •*- (Pq - 1) - {Rq +<5 0 + 1) 
m <- (Dq -Rq)—8q - 2 
mi=m 0 -(S 0 + 2) 

Thus the margin has decreased by 5o + 2. 

Similarly, at the end of a cycle, the divisor is one less than it was and the 
quotient has been increased by unity. Thus the difference at the start of the 
next cycle of the algorithm is: 

<5«- (Co + 1) - (Dq - 1) 

S^(Q 0 -D 0 ) + 2 

$ I = 0() + 2) 

At the end of the second cycle of the algorithm we have: 

m 2 = m | —(<51+2) 

'»2 =» i O-( (5 0 + 2 )-( <5 l + 2 ) 

m 2 = > n Q - ((5 q + 2) - (<5 q + 4) 
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In general, for positive values of n, we have: 

n 

m n = niQ - uSq - 2 ^ i 

i=l 

m n = ni() -uSq - 21 ~^n(n + 1) J 

mn = mQ—nSo>»nr -n 

m n = n?o “ n (^0 + 1) - n~ 

The margin between the remainder and the divisor gives some indication of 
how many more steps are required until there is the possibility of the 
remainder becoming equal to zero. This change occurs in line with the 
square of the number of steps of the algorithm. In particular, m„ can be 
equated to zero and the above quadratic solved for n. 

When the remainder is comparable in size to the divisor, the next step will 
either cause the remainder and the divisor to be equal, or the remainder will 
become greater than the divisor. If the remainder and divisor become equal, 
then the solution has been identified. If the remainder exceeds the divisor, 
the divisor is subtracted from the remainder once more. 

Recall that both divisor and quotient must be odd at a factorisation. The 
status of divisor and quotient are known in advance of applying the 
quadratic formula. Consequently, whether n is even or odd indicates whether 
or not the convergence might be a factorisation as shown in Table 1. 


Divisor 

Quotient 

n 

Possible 

Factorisation 

ODD 

ODD 

ODD 

NO 

ODD 

ODD 

EVEN 

NO 

ODD 

EVEN 

ODD 

YES 

ODD 

EVEN 

EVEN 

NO 

EVEN 

ODD 

ODD 

NO 

EVEN 

ODD 

EVEN 

YES 

EVEN 

EVEN 

ODD 

NO 

EVEN 

EVEN 

EVEN 

NO 

The prior state of the divisor and quotient together with the 

outcome of the quadratic formula determine whether the 

convergence of remainder and divisor offers any prospect of a 
factorisation. 

Table 1. 
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The behaviour of the remainder in the initial stages of the decrementing 
algorithm displays a characteristically quadratic shape as shown in Figure 6 
below. The data is taken from that shown in Figure 5 and is expressed in the 
form of the margin rather than remainder. 



200 195 190 185 180 


Divisor 

When the divisor and quotient are approximately equal, the margin decreases quadratically. 
The behaviour is less clear once the margin has reached zero for the first time. 
Figure 6. 


Once the margin has reached zero, the behaviour apparently changes. 
Comparison with Figure 5, suggests that the quadratic relationship may not 
hold beyond the first zero. There is, however, no restriction on the value of 
m n , and the value of n can be permitted to grow as required. Similarly, m 0 , 
and 5o are not restricted and can be based on the results of any choice of 
divisor. Consequently, new values of m 0 , and 5o can be established on the 
basis of the initial values and the outcome of the quadratic formula. A fresh 
quadratic may be formed and solved to find the next possible convergence 
region. 

Example 3 

In Example 2, the margin between the initial remainder and the divisor was 
44583 - 10434 = 34149. The initial difference between the divisor and the 
quotient was zero. Consequently at the start of the algorithm, the margin 
decreased by 2. Clearly there was no prospect of the remainder approaching 
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zero for a considerable number of repeats of the algorithm loop. Applying 
the quadratic formula to solve for n gives: 

-(go-Oo+l)± I(,Qo-D 0 +1) 2 +4.(D 0 -R) 


With the values above, this gives: 

-1+^1+4(34149) 


So in this instance, n = (-1 + 369. 2)12 = 184.3 Taking the nearest whole 
number for n, the divisor becomes 44583 - 184 = 44399. The quotient is 
incremented by the same amount, making it 44583 + 184 = 44767. The 
product of these two numbers gives 1987610033 which leaves a remainder 
of 44290. From Table 1, there was no prospect of a factorisation at this 
stage irrespective of the value of n since both divisor and quotient were 
initially odd. 

This is the first possible convergence of a pair of numbers to a remainder of 
zero and it is clear that with the 185th step, the remainder exceeds the value 
of the divisor. An extra single subtraction of the divisor with an associated 
increment of the quotient is then required to restore the correct remainder. 
Thus the Divisor becomes 44398, the Quotient = 44769 and the remainder 
becomes 44659 - 44398 = 261. Now the difference between the quotient and 
divisor has become 371, and the margin is 44137. 

Whereas n is the dominant term at the outset, the difference 5 between Q 0 
and D 0 soon assumes a greater significance. Implementation of this 
algorithm showed the expected rapid decline in the separation of the 
possible factorisations. 

As with the other algorithms, the introduction of the multiplier li restored 
performance whenever li was incremented. The calculations involved in the 
quadratic formula renders its use inappropriate whenever the values of n 
being returned are small. As an example, the first 300 iterations on the 
number 198762934254442223 was equivalent to some 364000 single 
decrements, an improvement of about 3 orders of magnitude. The first step 
was 18727. At the end of the 300 iterations, the step size was still 611, 
showing that the quadratic formula could continue to offer some benefit for 
many more iterations. 
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Considerations in the use of the algorithm 


Although many direct divisions or decrements are replaced by the quadratic 
approach, computers are sufficiently fast for the savings achieved by 
application of the quadratic algorithm above to be small even for numbers in 
the region of 10 18 . For much larger numbers, the benefits would be 
correspondingly greater. 

The step size calculated by the quadratic formula is rounded down to the 
nearest integer. As the step size diminishes, the algorithm loses its potency. 
The algorithm comes to a complete halt when the step size is rounded down 
to zero. It is therefore necessary to alternate between the quadratic algorithm 
and one of the other approaches. 

The Newton Raphson Method 

Calculations involved in the standard quadratic formula are somewhat 
onerous, but it is worth noting that the curve of the step size is smooth. The 
Newton Raphson method can be applied using the step size established 
between one pair of potential roots as the starting point for the solution of 
the next quadratic. The required accuracy is to the nearest integer, so the 
Newton Raphson method could be expected to converge rapidly. The first 
iteration would be given by: 

m r D r R i 

s j =Q j~ D j 

fi n j) = mj - n j(dj + 1) - nj 
f(n j ) = -2n j -(S j +l) 

A n j) 

nj+l=nj ~7^) 

The Newton Raphson formulation avoids the need for evaluating the square 
root, and the multiplications to be carried out are more modest in size than 
in the quadratic formula. The result of the division which is carried out 
raises a potential problem. In general the result of this division will not be 
an integer and it would be appropriate to proceed to the second and 
subsequent iterations using floating point methods. Additionally, the test for 
convergence requires a redundant application of the algorithm. 
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Normally, this is of no consequence since the Newton Raphson method 
would usually be applied to obtain a root to a fairly high accuracy. In such a 
scenario some number of iterations would be expected, and the additional 
loop required to confirm convergence would be little consequence 
computationally. 

In the present application, the initial estimates are expected to be fairly 
accurate. The division could therefore be rounded to the nearest integer. 
This would avoid the need for floating point arithmetic, but there is an 
attendant risk that the iteration would then become oscillatory. It should, 
nevertheless be possible to reach a close estimate for the value of n j +\with a 
single application of the Newton Raphson method. Thereafter the 
coefficients could be adjusted using the decrementing algorithm. 

Finite Differences 

An alternative approach alleviates the problem of root finding, to an extent, 
through judicious application of finite differences. The smooth nature of the 
curve of separation of convergences, permits the use of a simple 
extrapolation algorithm based on finite differences to predict the locality of 
the next potential zero. Some actual step sizes for factorisation, together 
with their first and second differences are illustrated in Figure 7. 



Note that although successive convergence intervals remain large, the first and second 
differences rapidly converge towards zero. 

Figure 7. 
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Table 2 contains some of the data used to construct Figure 7. 


Convergence First Second 

teration Interval Difference Difference 

1 

18727 



2 

9494 

-9233 


3 

7022 

-2472 

6761 

4 

5839 

-1183 

1289 

5 

5107 

-732 

451 

6 

4595 

-512 

220 

7 

4212 

-383 

129 

8 

3912 

-300 

83 

9 

3668 

-244 

56 

10 

3464 

-204 

40 

11 

3291 

-173 

31 

12 

3141 

-150 

23 

13 

3011 

-130 

20 

14 

2896 

-115 

15 

15 

2792 

-104 

11 

16 

2699 

-93 

11 

17 

2615 

-84 

9 

18 

2538 

-77 

7 

19 

2467 

-71 

6 

The intervals between potential factorisation points decline in a manner which 
can be predicted with reasonable accuracy using finite differences. 

Table 2. 


Example 4 

By reference to table 2 it can be seen that after 17 iterations, the step size is 
2615. The first difference was -84 and the second difference had fallen to 9. 
This means that an extrapolation based on the previous step size and the first 
difference alone would be within 10 of the next potential factorisation. 
Using both differences in this example, the next step size is given by 
2615-84 + 9 = 2540. The value determined by the quadratic formula was 
2538. Simple arithmetic operations have produced a result which is within 
0.1% of the quadratic formula at a fraction of the computational cost. This 
estimate would subsequently be adjusted using the decrementing algorithm 
to locate the exact point of closest approach and identify whether it is a 
factorisation. The actual step size is then updated along with the most recent 
first and second differences prior to estimating the next step. 
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Although it is possible to run the decrementing algorithm in reverse, it 
seems simpler to employ the first difference on its own. This gives a 
consistently low estimate for the separation of convergences and the 
decrementing algorithm can then be applied without concern over the 
possibility of missing a factorisation. 

Another method for estimating the next convergence interval is to divide the 
divisor D by the difference 5. With every unit decrement in the divisor D, 
the margin between D and 5 decreases by an amount which is at least as 
large as 5. When 5 is very small, the error is gross and the quadratic 
approach is to be preferred. As 5 increases in size, the accuracy of the 
estimate improves dramatically. The estimate is always high since no 
account is taken of the progressive decrease in the size of the divisor and the 
corresponding increase in the size of the quotient which occurs in the 
intervening steps. In the above example, by the time D/d ~ 1000 there was an 
error of about 1%. The estimate may then be reduced by a small amount 
followed by a switch to the decrementing algorithm until the nature of the 
point of convergence has been determined. Thereafter a fresh estimate of the 
number of decrements to the next convergence can be made. 

Returning to the Decrementing Process 

Having determined the step size and evaluated a new estimate for the divisor 
and quotient, the corresponding remainder needs to be determined by 
calculation. Simplistically, this involves either a division of the target 
number by the new divisor, or multiplication of the new divisor and 
quotient. 

The problem can be reduced slightly by noting the relationship which the 
new divisor and quotient have to the previous values employed. Prior to 
carrying out adjustments, the divisor and quotient are both changed by the 
value of n. Let: 


D r D i- 1 -” 

Qj = Qj -1 +n 
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Then: 


R j = T ~ (D j~ 1 _n X2/-l +”) 

Rj = T-Dj_ x Qj-\ +n(Qj_\ - Dj_ { ) + n 2 
Rj=Rj-i +n(Qj _i +«) 


In this formulation, the product DjQj has been replaced by a smaller 
multiplication. 


Algorithm Flow 


It is evident that although the quadratic approach is beneficial at certain 
stages of the factorisation search. The finite difference estimate or the direct 
division estimate can then serve as a surrogate for the quadratic formula. 
However, these approximations also incur a certain amount of overhead. It 
is clear that as the step size diminishes, there comes a point at which it is 
more efficient to persist with the decrementing process. This situation will 
persist until the divisor has diminished to a point at which li requires to be 
incremented. Following the increment in p a return to the quadratic formula 
might be warranted. 
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A Third Quadratic Search Algorithm 


Observe that the quadratic formula and its alternatives only seek to identify 
possible regions of convergence between the values of the divisor and of 5. 
Whether or not a convergence is a factorisation has to be determined 
explicitly. As the convergences become more closely spaced, the question of 
whether or not any of them provide a zero remainder is worthy of further 
examination. In this context, the relationship between the divisor D and the 
interval between convergences, n becomes important. It was noted earlier 
that in general, there is no requirement for the margin to be set equal to zero 
in order to determine convergence intervals. In particular, the margin can be 
permitted to become negative. When a factorisation is found, the margin 
will be an integer multiple of the value of the divisor at that point. So, for 
some leN the quadratic relationship becomes: 

m n = -k(DQ - n ) = m q - h(Sq + 1) - n 2 
- kDQ + kn = m o - n(S q + 1) — n~ 

.'. rr + (Jq + k + 1 )n - (kDQ + hiq) = 0 

Applying the quadratic formula in this case gives: 

— (<5q+At+ 1 )± J (S Q+k+ \ +4 .(kD Q+m q) 


Each of the roots, when subtracted from D n , give a potential factor of the 
target number. The only permissible values of k are those which give rise to 
a discriminant which is a perfect square. Expanding the discriminant leads 
to: 


(Sq + l) 2 + 2k(S 0 + 1) + k 2 + 4D 0 k + 4m 0 
.'. k 2 + 2(d 0 + 1 + 2 D 0 )k + ((Sq + l) 2 + 4m 0 ) = L 2 
Where L e N. The above may be expressed in the form: 

(k + y) 2 = L 2 - <f) where: 

y = (Sq + 1 +2Dq) and: 

0 = (Sq + l) 2 +4 mQ-y 2 
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The above equations set upper and lower bounds on L since k + y is real and 
0 <k<T. Consequently the range for L is given by 
Jy 2 + (j> <L< J(T+y) 2 + </>. Only those values of L which give rise to a 
perfect square are permissible. 

The analysis is essentially equivalent to the original quadratic search 
algorithm. The value of k may be incremented, giving rise to a possible 
value of L 2 which is tested to establish whether or not it is a perfect square. 

Increments may be made to k and L alternately in the manner of the r|5 
algorithm. Thus k may be incremented causing the LHS to be greater than or 
equal to the RHS of the above equation. If it is greater, then L is 
incremented until the RHS is greater than or equal to the LHS. Equality ends 
the search. 
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A Multiplication Approach to Factorisation 

In a long division process, each step is actually a test multiply followed by a 
subtraction and the introduction of the next digit of the target number. The 
quotient and remainder can only be known at the end of this process. Each 
new test factor requires a fresh long division. The possibility of using test 
factors to multiply together rather than dividing the target number by a 
single test factor is worthy of consideration. 

Normal Multiplication 

When two numbers are multiplied together, it is customary to follow a 
particular sequence of steps. The first number is multiplied throughout by 
the first digit of the second number. Multiplication of the first number is 
then repeated using the second digit of the second number. This process 
continues using each of the digits of the second number in turn. Finally, the 
results of all the multiplications are added together to obtain the required 
result. In a normal situation, the calculation is carried to completion before 
any operations are carried out on the result. As with division, multiplying 
large numbers together requires a considerable number of steps and the time 
taken is commensurately large. 

A Modified Multiplication Approach 

During a multiplication process, the digits in the product can be made 
available starting with the least significant digit. Comparison between the 
product and the target number can proceed on a digit by digit basis, starting 
with the least significant digits. 

Suppose a pair of test factors T1 and T2, are to be multiplied together in an 
attempt to regenerate a target number. The first step is to multiply together 
the least significant digits of the two test numbers. If the least significant 
digit of the result is not the same as the least significant digit of the target 
number, then the test fails and a different pair of test numbers can be 
selected. If the first digit is correct, then the test continues. The 2nd least 
significant digit of T1 is multiplied by the least significant digit of T2, and 
the least significant digit of T1 by the 2nd least significant digit of T2. The 
least significant digits of the results are added together and to any carried 
digit from the first test. The least significant digit of the result is compared 
to the 2nd least significant digit of the target number. If they are different, 
the test fails. Otherwise, the test proceeds as above to the third digit and so 
on. Evidently, most attempted multiplications will fail long before the full 
multiplication of the test numbers has been carried out. As soon as the test 
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fails, new test numbers can be substituted and the process repeated. The 
procedure is illustrated in Figure 8 below. 



Test Fails at 
the third step 



Progressive multiply and compare procedure. 

Figure 8. 


As an alternative to the above digit by digit multiplication, some of the least 
significant digits of both numbers can be fully multiplied together. 
Comparison is then made between the lowest order digits of this product and 
those of the target number. 

The number of digits extracted from the product and used for comparison 
purposes has to be carefully limited. Higher order digits are influenced by 
the more significant digits of the test factors. Figure 9 illustrates the 
approach. 
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3 12 1 (£T)3 V 


1 9 4 


digits of result 


y 

Higher order digits are 
unsuitable for comparison 


Comparison can be made between the target number and the lower order digits from the 
result of a multiplication which uses only the lower order digits of the test factors. 

Figure 9. 

Either of the above arrangements may be implemented using standard 
programming facilities. Although it could not be expected to be as efficient 
as a specially tailored algorithm, the procedures offer potential for 
performance improvement when compared to a full multiplication prior to 
carrying out a comparison between the result and the target number. 
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Reverse Division 


The multiplication process described above suggests a potential means of 
factorisation. Consider the least significant digit of the target number. Only 
certain numbers multiplied together can give rise to any particular digit as 
illustrated in Table 3. 


Possible Product Pairs 

i 

lxl 

3x7 

9x9 


2 

1x2 

2x6 

3x4 

4x8 | 6x7 | 8x9 | 

3 

1x3 

7x9 


4 

1x4 

2x2 

2x7 

3x8 

6x4 

6x9 | 8x8 | 

5 

1x5 

3x5 

5x5 

7x5 

9x5 


6 

1x6 

2x3 

2x8 

4x4 

4x9 

6x6 | 7x8 | 

7 

1x7 

3x9 


8 

1x8 

2x4 

3x6 

4x7 | 6x8 | 2x9 | 

9 

1x9 

3x3 

7x7 


0 

1x0 

2x0 

3x0 

4x0 5x0 6x0 7x0 8x0 9x0 0x0 

The least significant digit of a product, can only arise from 
certain combinations of numbers. 

Table 3. 


The simplest application of Table 3 is in establishing whether or not a 
number is a perfect square. The usual subroutines available in programming 
languages calculate the square root of a number. When writing a program to 
determine if a number is a perfect square, a straightforward approach is to 
calculate the square root and then test the result to establish whether or not it 
is an integer. In reality, this incorporates a lot of unnecessary calculation. 

Instead of checking to determine whether or not the square root is an 
integer, the problem can be posed in a slightly different way. 

/ \ 

Is there some integer which when 
multiplied by itself will result in a 

product equal to the number under 

V___ J 

Examination of Table 3 shows that only certain least significant digits can 
be generated by multiplying two identical numbers together. It has already 
been noted that an odd perfect square can only result from an odd square 
root, an even perfect square can only result from an even square root, and no 
perfect square can end in 2, 3, 7, or 8. 
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Example 5. 


Assume that the possible perfect square to be tested is 13826769. Start with 
the least significant digit. Examining row 9 of Table 3 shows that if the 
square root is an integer, then it will end in either a 3 or a 7. Figure 10 
examines the first of these options. 


PoSSi UI " t 1 -.-*’'..- 1 c... _ /A A\n O /tA 1 7 c. I I : 



V 6 9 

The possible perfect square must be the result of the 
multiplication of a pair of identical integers. 

Figure 10. 


The second least significant digit of the possible perfect square is considered 
next. There is no carry digit from multiplying 3x3, and the product of the 
possible square root with itself must produce a 6 in the second column to 
match the number under test. Thus the second least significant digit of the 
possible square root can only be 1. This choice gives rise to a carry of 1 into 
the third column of the product of the possible square root with itself. 

Progressing to the third digit of the possible perfect square, the square root 
requires to produce a 7 in this column. There is a carry of 1, so when taken 
as a product with the least significant digit of the possible square root the 
number chosen has to produce a further sum of 7 - 1 = 6. In other words, the 
product with the least significant digit has to be 6/2 = 3. Thus the third digit 
of the possible square root can only be 1. This choice provides the correct 
third digit and results in a carry of 2 into the fourth column. 

The fourth digit of the possible perfect square is 6. This requires that the 
fourth digit, x of the possible square root be such that 2(3x) + 2 = *6. 
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The asterisk represents a possible carry digit. Referring to Table 3, 
3x=*2=>x = 4. Although this produces the correct result for the fourth 
digit of the possible perfect square, the product of the possible square root 
has carried to 8 digits with a value which is greater that the possible perfect 
square. The process has therefore failed. 

Having established that the possible square root does not end in 3, the 
search is repeated for a number which might end in 7. In this case, there is a 
carry of 4 into the second column. The number chosen for the second 
column of the possible square root then has to satisfy 2(7x)+4 = *6. The 
choice from Table 3 must then be such that 7x=*l. The only possible 
number is 3, which gives a carry into the third column of 4. The possible 
perfect square has a third least significant digit of 7. Therefore the choice of 
third digit for the possible square root is defined by: 

2(7x) + 4 = *7 

2(7x) = *3 

There is no value of x which can satisfy this relationship since the right hand 
side ends in an odd number. The process again fails. All the possible routes 
to an integer square root have been exhausted. The conclusion is that the 
original number is not a perfect square. 

In the case of square roots, the digits of the quotient and divisor are identical 
which simplifies the problem. 

Table 4 is an alternative representation of the information in Table 3, which 
can be used to find a quotient given a target number and candidate divisor. 



i 

2 

3 

4 

5 

6 

7 

8 

9 

0 

1 

i 

- 

7 

- 

- 

- 

3 

- 

9 

- 

2 

2 

1;6 

4 

3;8 

- 

2;7 

6 

4;9 

8 

- 

3 

3 

- 

1 

- 

- 

- 

9 

- 

7 

- 

4 

4 

2;7 

8 

1;6 

- 

4;9 

2 

3;8 

6 

- 

5 

5 

- 

5 

- 

1;3;5;7; 

9 

- 

5 

- 

5 

- 

6 

6 

3;8 

2 

4:9 

- 

1;6 

8 

2:7 

4 

- 

7 

7 

- 

9 

- 

- 

- 

1 

- 

3 

- 

8 

8 

4;9 

6 

2;7 

- 

3;8 

4 

1;6 

2 

- 

9 

9 

- 

3 

- 

- 

- 

7 

- 

1 

- 

0 

0 

0;5 

0 

0;5 

0 

0;5 

0 

0:5 

0 

0 

Using a possible divisor for a known product, a corresponding quotient may 
be constructed by successive application of the reverse division table. 

Table 4. 
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Starting with the target number the table is used as follows. 

The least significant digit of the target number is located on the left column, 
and the least significant digit of the divisor along the top row. The entries 
where the row and column intersect are the possible least significant digits 
of the quotient. If there is no entry, then the potential divisor will not divide 
the target number exactly. Thus, for example, no number ending in 7 can be 
divided by any number ending in 8. Similarly, if the target number ends in 4 
and the potential divisor ends in 6, there may be an integer quotient, and if 
so it will end either with a 4 or with a 9. 

Having identified the least significant digits for any possible quotient, the 
penultimate digits for these quotients can be determined. This is done by 
taking the least significant digit from the first of the possible quotients 
determined as above and multiplying by the potential divisor. The product 
gives the appropriate carry value which is subtracted from the penultimate 
digit of the target number. Applying the division table to the result 
determines possible penultimate digits for the first possible quotient. This 
procedure is repeated for each possible quotient in turn. 

In the case of an even divisor or one ending in 5, the process fails if a step 
gives rise to an entry for which no product is possible. For divisors ending 
in 1, 3, 7, or 9, failure will only occur with overflow at the end of the reverse 
division process. It can be seen that the technique creates an increasing list 
of possible quotients which eventually either fail or overflow. 

Reverse division is not restricted to determining square roots or finding 
quotients given test divisors. It is possible to obtain a factorisation directly 
by application of the principle. 

In the general case of the quotient and divisor being distinct and it is 
necessary to consider each of the pairs of numbers which could have given 
rise to each digit in turn. As each more significant digit of the target number 
is considered, the number of factor pair candidates increases. 

Since all divisors of interest will be odd, failure will only occur at overflow. 
Thus the method would not be practical for complete factorisation. 
However, it could serve to generate a list of acceptable ending groups of 
numbers, which could be used in conjunction with prime number generation 
to eliminate candidate divisors. 
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Consideration of the discriminant and the relationship 
between pand 8. 

For any quadratic equation, the nature of the roots is known from the value 
of the discriminant: 

+4/iT 

The factorisation problem as formulated herein, ensures that the 
discriminant is always positive. This guarantees that the roots are not 
complex. When the discriminant is not a perfect square, the roots are 
irrational and the solution is of no interest. For any values of 8 and p giving 
rise to a discriminant which is a perfect square, the roots are rational 
numbers. 

As already noted, for any monic polynomial with integer coefficients, the 
only rational roots which can exist are integers. Furthermore, such roots are 
factors of the constant term. In the basic form of the first quadratic 
algorithm, p is required to be unity and the discriminant is always required 
to be a perfect square. The term 5 2 is also required to be a perfect square to 
ensure that the only coefficients considered are integers. This guarantees 
that any roots found are integers and therefore factors of the target number. 

Thus, although the formulation is constructed in terms of the divisor D, the 
second root returned turns out to be the associated quotient Q. 

In the improved algorithm, p is not required to be unity, but is nevertheless 
an integer. The polynomial is no longer monic. Consequently there is no 
theorem for the general case which guarantees that any real roots will be 
integers. However, the condition that the discriminant be a perfect square 
guarantees that the roots are at least rational. Insisting upon the term 5 2 
being a perfect square also guarantees that all the coefficients are integers. 

The relationship between the roots and the coefficients of a quadratic 
equation is well known. Denoting the roots by a and p we have: 

(a + fi) = -jj 

“P=jr 

From the product of the roots above, it is certainly possible that there may 
be integer roots a and p when p divides T exactly. However, in that case, p 
is already a factor of T and the problem has been reduced. 
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Should p divide T exactly, and also divide 5 exactly, a and P can both be 
integers because the quadratic has been reduced to a monic polynomial with 
integer coefficients by dividing throughout by p. 

On the other hand, it can be seen that if p does not divide T exactly, then 
there can be no solution with both a and p integers. Nevertheless, it is 
possible for one of the roots to be an integer. 

From data produced by the program, it appears that when the discriminant is 
a perfect square, the positive root a is always an integer and the negative 
root p is always rational. From the product of roots given above, it can be 
seen that the relation between the Divisor D, quotient Q and the roots is 
then: 


D = a 


Q=t*P 


In the factorisation problem both the discriminant and 5 are required to be 
integers. The roots a and p therefore assume the form: 

_ ~^+ ^ 2 +4/<r _ ^_ § 

a 2ju 2jii 

-5- jd 2 +4 f .iT _(2 + <S) 

P ~ 2n ~ 2 ju 

Since A. is always greater than 5, a is always positive and p negative. The 
roots are clearly rational. However, a will not be an integer unless X-3 is a 
multiple of 2p. Similarly for p with respect to X + S. Applying the second 
relationship between quadratic roots and the equation coefficients, we have: 

(A-w5)(A-(5) _ 

¥ - 1 

Since T is odd, it is required that the roots a and p when multiplied together 
produce an odd integer. Thus X-S and X+8 must each be odd with the 
product being an odd multiple of 4p. 

This in turn means that whenever 8 is odd, the discriminant must be even, 
and when 8 is even, the discriminant must be odd. 

Recall that the quotient is defined in terms of the divisor thus: 

Q={iD+S 
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The value of p has a relationship to 5 to the extent that even values of p 
exclude even values of 8 in the factorisations. This is to be expected. Since 
the target number is odd, both factors must be odd irrespective of whether 
they are composite or prime. If p is an even number, then the product of p 
taken with D will also be even. Thus Q can only be odd provided an odd 
value of 5 is added to this product. Conversely, when p is odd, its product 
taken with the D will also be odd. Consequently, 5 is required to be even to 
ensure that Q will be odd. 

When applied to the quadratic search algorithm, the above relationships 
provides extra tests for 5 2 . It has already been noted that no perfect square 
can end in 2, 3, 7, or 8. For even values of p which occur during operation 
of the algorithm, 5 is required to be odd. When 5 is odd, so too is 5 2 . 
Consequently it can only end in 1, 5, or 9. For odd values of p, 5 2 is required 
to end in 0, 4, or 6. 

The test can be implemented by counters modulo 10 with suitable offsets for 
1, 5 and 9 used when p is even and with offsets for 0, 4, and 6 in those cases 
where p is odd. As an alternative to using modulo 10 counters, there is a 
cyclic pattern in the values of the final digit of 5 2 which might be exploited 
to allow an appropriate step size in X. It would then be possible to avoid 
generating those values of 5 2 which would be rejected by the modulo 10 
tests. 

In the first quadratic algorithm for the factorisation problem, the 
discriminant assumes a particularly simple form. As noted, in the early 
stages of incrementing X, when 5 is small, many values of 5 are skipped 
over. It was then shown that the search method using q and X obviated the 
need for explicit testing of the value of 8 2 . 

The second quadratic algorithm permits considerable savings in the range of 
values where 5 is small. By employing the finite difference approximation 
for the quadratic formula, this range was extended. 

Eventually, however, recourse must be made to some other method since the 
second quadratic algorithm eventually stalls. Both the decrementing 
algorithm and direct division have been considered as alternatives after the 
transition. The i], X approach to the first quadratic algorithm can also be 
employed for this phase as follows. 

The values of 5 and p at the point of transition can be determined. From 
these, and the value of the target number T, an initial value of X can be 
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calculated. The value of q is made equal to 5. Since 5 will be fairly large at 
the transition, conditions are appropriate for using the r| and X approach. 

An upper limit on the value of i] which is required prior to a return to the 
second quadratic algorithm can be derived as follows. Recall that an 
increment in p occurs when: 



At this point, 5 is reset to zero and Q = (ju + 1 )D. When 5 becomes equal to 
the above value of D an increment in p is required. Thus the q and X search 
may be applied between the points at which the marginal quadratic search 
stalls and the next increment of p. Since even values of p require odd values 
of 8 and vice versa, q may be incremented appropriately. 

Since D must always be an odd integer, it can be seen that k-8 must always 
be an odd multiple of 2p. It is also clear that when 5 is odd, so too is X and 
when 5 is even, X must be even. Consequently, in the q and X search, as 
each new value of q is tested, the next value of X is determined by the value 
of 2p. So in fact, although the range which X covers is large, as the value of 
p increases, an increasing number of the values of X are omitted from the 
search. Furthermore the higher values of p coincide with smaller ranges in 5. 

It has also been noted that elimination of the possibility of small prime 
factors in T would be carried out prior to any of the search procedures 
described. Consequently, since Q=juD+d , then at least for the smaller 
values of p, those values of 5 which are divisible by p or its factors (other 
than 1) are also immediately excluded from a possible factorisation. 
Otherwise, Q and hence T would be divisible by p or one of its factors. This 
possibility has already been eliminated by the initial tests. 
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Tables of Prime Numbers and the Sieve of Eratosthenes 

Despite the speed and simplicity of the decrementing algorithms described 
here, every odd integer assumes the role of divisor in the process. Most of 
these are composite numbers, so it would seem preferable to skip to the next 
prime divisor. Creating a lookup table of prime numbers is a straightforward 
matter. The step size C between successive primes is readily derived from 
such a table and may be saved as a separate table. Having established a 
value of divisor which is a prime, the corresponding quotient and remainder 
can be determined. In order to arrive at values for the quotient and 
remainder corresponding to the next value of divisor, the following 
procedure may be applied: 

D ■*- D-C 
R+-R + CQ 
Q *- Q + R div D 
R <- R mod D 


Example 6. 

Target Number T = 2631943 
Initial Divisor D = 307 
Initial Quotient = 8573 
Initial Remainder = 32 

Assumed step size C = 24 
New Divisor Dj = 307 - 24 
D, = 283 

Interim Remainder R’ = 32 + ( 24 x 8573) 

R' = 205784 

New Quotient Qj = 8573 + (205784 div 283) 

Qj = 8573 + 727 
Qj = 9300 

New Remainder Rj = 205784 mod 283 
R, =43 

The multiplication and division required in these calculations are slightly 
smaller than would be required to divide into the target number. 
Consequently they will be faster in execution. The approach bridges the gap 
between a direct division approach and the decrementing techniques. 
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A lookup table eliminates all composite numbers in the range which it 
covers. Beyond this range, a sieve of Eratosthenes can be constructed to 
eliminate a high proportion of all composite numbers. 

A Sieve consists of modulo P counters, one for each prime P in the sieve. 
The counters operate synchronously. Each of the sieve primes in turn is 
divided into the initial divisor and the associated modulo P counter 
initialised with the remainder. Once initialisation is complete, the divisor 
register and sieve counters are decremented synchronously. 

The output from any sieve counter will equal zero when the divisor assumes 
a value which is divisible by the prime number associated with that counter. 
Consequently, if any one or more sieve counters = 0, the number in the 
divisor register is definitely composite. If all the sieve counters ^ 0 then the 
number in the divisor register is prime to all the numbers in the sieve and 
may be a prime. 

The software implementation of a sieve is cumbersome. However, a set of 
wired OR modulo N counters is a practical alternative, and could readily be 
constructed for a few hundred prime numbers. 
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The Nature of Composite Numbers 


It has been noted that in the algorithms described, many of the numbers 
serving as divisors are known to be composite. Ideally, divisors should be 
prime numbers. It is worth considering the effect of using divisors which are 
composite numbers. 

In the case of a highly composite target number, there are several possible 
divisors. For example, consider a number being tested which is the product 
of 6 primes. At most one of the primes is larger than the square root of the 
number. If that is the case, one factor is the product of all the primes less 
than the square root. Thus, in this example, if one factor is greater than the 
square root, then there are 5 prime numbers, all of them less than the square 
root, which will identify a factorisation of the original number. However, 
there are also 5 C? + 5 C 3 + 5 C 4 + 5 C 5 = 26 composite numbers which will lead 
to a factorisation. These composite numbers are distributed over the range 
from unity up to the square root of the target number. Thus irrespective of 
where in the range, up to the square root of the target number, that a search 
algorithm was started, there would be at least some prospect of finding a 
factorisation. Since the nature of the target number is not generally known in 
advance, incorporating composite numbers as divisors may not be such a 
serious drawback as might at first be assumed. 

Conversely, in a composite number T which consists of the product of n 
primes, at least n - 1 of these primes are smaller than ff . The largest factor 
smaller than ff will then be the product of m<n-\ primes. 

Consider the said largest factor as a new target number Ti. By application of 
the same reasoning, it follows that at least m -1 of the primes which 
comprised the original target number are smaller than fff and that at most 
one of the primes p : JT\ <p<T\. Thus, at most two of the primes from the 
initial target number are greater than iff. 

Taking the argument further, consider the case in which all the factors are 
equal to iff . The product of the factors equals T. Suppose now that any 
single factor in this set of factors was greater than iff . Then the product of 
all the factors would exceed the value of T. It follows therefore, that for any 
composite number T consisting of n factors, there exists at least one factor 

Pi<fT. 

In the example above, where the target number consisted of 6 primes, one of 
the prime factors must be less than iff . 
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It is evident that this characteristic of composite numbers favours an 
upwards search as a means of finding a factorisation. 

The relationship between the number of prime factors which comprise the 
target number and the corresponding root of the target number is an 
important characteristic. It offers a progressive test for determining the 
maximum number of prime factors which can be contained in a given target 
number. For example, assuming that an upward search has failed to identify 
any factor less than !fr. Then the target number can have at most k - I 
prime factors. 
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Diophantine Equations 

Equations which require solutions consisting of sets of integers or rational 
numbers are classed as Diophantine equations. 

Pythagoras’ equation can be applied to solve any right angled triangle, 
irrespective of whether the sides have lengths which are integer in length. 
Clearly by choice of a suitable unit of measurement, at least one of the sides 
can be represented by an integral value. However, restricting the problem to 
a search for triangles which have integral values for all three sides brings it 
into the class of Diophantine equations. 

The Pythagorean numbers are integers which satisfy the equation 
x“ +y“ = V~. A partial solution is given by the substitution: 

x = 2n + 1 

y = 2«“ + 2 n 

z = 2+ 2n + 1 

A general solution to the problem of finding integer triples which satisfy the 
Pythagorean equation is given by the Pythagorean numbers defined by: 

2 2 
x = m - n~ 

y = 2 mn 

2 2 
z = m + n 

The general solution offers flexibility in the production of Pythagorean 
numbers. On the other hand, given values of x, y and z, the identification of 
the values of m and n is not so convenient. The first set of equations do not 
give rise to as complete a range of triples as does the second set. However, 
there is only one parameter instead of two. Thus given any odd value of x, 
the value of n can be determined. From that the values of y and z follow 
immediately. 

In the factorisation problem we have that X 2 = 5 2 + 4T, and it is natural to 
think of X as the length of a hypotenuse and the other numbers as the lengths 
of the other sides of some right angled triangle. The relationship between 
these numbers is reminiscent of the 3, 4, 5 triangle or the 5, 12, 13 triangle. 
It is clear that setting n = tin the first group of equations above gives rise to 
the 3, 4, 5 triangle while n = 2 returns the values of the 5, 12, 13 triangle. 
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The comparison is not exact since 4T is not a perfect square for factorisation 
problems of interest. The corresponding length is 2 JT which is therefore an 
irrational number. The actual equation of interest is: 


This form matches right angled triangles of the type 1,2,73". More 
specifically, the similar triangle 2, 4, 2 73" matches the form of interest. 

For the factorisation problem, x = 4 T and is clearly even. From the 
requirement that both of the factors be odd it follows that y s <5 and - = /. are 
both even numbers. 

Evidently, any pair of even numbers can serve for y and z and will give rise 
to some value of x viz: 


Let : y = 2a\z = 2b 
Then : x = Ab~ -Aa~ = A(b~ -a 2 ) 

.-. T= b 1 - a 1 

While the generation of T is trivial, the inverse problem of finding values of 
y and z which generate a given value of T remains. It can also be seen that 
since T is odd in the factorisation problem, one of the pair of numbers a and 
b defined above must be odd while the other must be even. So we may 
write: 


b = 2 n 


a = 2 « - 1 


So that z = 4 n and y = 4n - 2 

.'. x = 16n^ - (An — 2)^ = 16 - {16 - 16 n + 4} 
x = 16n - 4 = 4(4« - 1) 
or : T = An - 1 

The form satisfies the requirements of the equation and relies on the single 
parameter n. Consequently, given any one of the three terms, the other two 
can be determined. However, the above definition of x, y, and z is not 
particularly useful since the values returned lead to a factorisation with a 
divisor of unity and a quotient of T. 
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It has already been shown that with S = T- 1 the determinant +4T is 
always a perfect square and gives X = T+ 1. This guarantees that at least one 
right angled triangle with sides <5 , X , 2 JT can always be constructed for any 
value of T. The first few of these triangles are illustrated in Figure 11. 



I- 12.000 -1 

The first few triangles in which the hypotenuse and base 
are integer lengths. The third side is of a length equal to 
twice the square root of an odd number. 

Figure 11. 

The included angle is determined by the ratio of the sides: 

CosW) = j = j^ 

Successive even values of T give rise to the sequence: 

The sequence for even numbers is not of interest in the factorisation 
problem, since even values of the target number can always be divided by 2. 

Considering successive odd values of T, the possible values of 6 conform to 
the sequence: 

CosW)e{±; 

This may be seen by reference to Figure 10. The above sequence defines 
one set of permissible values of the ratio S/L 
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Taking as an example the case of T= 9, an additional solution exists in 
which <5 = 0,1 = 9 and for which the triangle collapses. This type of 
additional solution will exist for all numbers which are perfect squares. 

Numbers which can be factorised, also have additional solutions. Each 
additional solution results in the construction of a triangle in which the value 
of 0 is given by: 

Cos(9) = t = 

a (q + d ) 

Since : AT = + — <5~ = (q + d)^ - (q - d)^ 

AT = q- + 2 qd + d~ -q^ + 2qd - d^ 

T = qd 

Thus the quotient and divisor are related in a simple manner to the lengths 
of the sides of a right angled triangle with sides of length 2 JT , 8 and X 
This is illustrated in Figure 12. 



Right angle triangles corresponding to the three prime factors 
of 105. (3, 5, 7). The base length is (q - d) and the hypotenuse 
is (q + d). The third side is twice the square root of 105. 


Figure 12. 


Multiplying top and bottom by d gives: 


50 



Cos(O) = 


d(q—d ) T~d~ 
d(q+d) T+d 1 


(1 <d<jT) 


Thus given any value of T, it can be seen that for any factorisations which 
exist, the value of 6 belongs to the sequence: 


Cos(d) = 


(T-l) ' (7-9) . (T- 25) . 
(r+i) ; (t+ 9) ; (r+25)’’".’ 


All of the triangles formed are right angled, and from elementary geometry 
it can also be seen that the factorisations correspond to a series of circles as 
illustrated in Figure 13. The radius of any circle is half the length of the 
hypotenuse defined by the chosen factorisation and is equal in length to XU. 
All of the centres lie on the perpendicular bisector of the side, the length of 
which represents 2 JT. The distance between that side and the centre is 3/2. 



The illustrations given are in terms of integer solutions, but X and 5 are 
defined in terms of d and q. The only restriction on the values of d and q is 
that they satisfy T = qd. However, the existence of a partial solution defined 
in terms of a single parameter is encouraging. The possibility of a set of 
such solutions addressing a variety of factorisations is a tantalising prospect. 
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An Analogue Approach 

The difficulty in finding a solution to the factorisation problem lies in the 
fact that any real number can be divided into the target number and will 
return another real number. The divisor and quotient real number pairs may 
be irrational, rational or within the integer subset of the rationals, but 
nevertheless an infinite number of factorisations exist for any target number. 

The only factorisations of interest are those in which both factors are 
integers, but there is no obvious mathematical constraint which can be 
placed on the equation DQ = T to force the only possible solutions to be 
integers. 

Instead, the integer solution aspect has to be approached mechanistically, by 
ensuring that the only numbers tested are integers. This, in turn requires a 
stepwise approach, the complexity of which is inevitably related to the size 
of the target number. 

If there is to be an analytical solution which is essentially independent of the 
size of the target number, some means of imposing the integer constraint is 
required. 

In nature there is an extensive range of examples where integer solutions are 
the only ones possible. Quantum Mechanics relies heavily on integer 
solutions to describe optical, atomic, and nuclear phenomena. 

There are also much more commonplace examples such as musical 
instruments in which harmonics are related to the fundamental note by 
integer multiples. 

Factorisation by a Wavelength Analogue 

By definition, any divisor of the target number fits into the target number an 
exact number of times. Thus in a sense, a factor is equivalent to a harmonic. 
If the value of the target number is regarded as the wavelength of a 
fundamental, then a divisor D corresponds to the Qth harmonic of the 
fundamental, where DQ = T. 

In the case of a square wave, every odd harmonic is present and the wave 
may be described by its Fourier Series: 

00 

x(f)=X sin(fa»ot) ( kodd ) 

k= 1 

So a square wave signal provides a means of generating all the possible 
divisors for any number continuously and simultaneously. 
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Any square wave will posses the full range of harmonics irrespective of the 
fundamental frequency used. However, the fundamental frequency chosen, 
corresponds to some integer wavelength, and only certain of the harmonics 
can then have wavelengths which also correspond to integer values. 

Furthermore, although a perfect square wave contains an infinite number of 
harmonics, the highest harmonic required is determined by the smallest 
integer divisor to be tested. 

In principle it is possible to construct a set of resonators, each of which is 
tuned to an integral wavelength. This set of resonators would then be driven 
by the square wave. If any harmonic had an integral wavelength equal to one 
of the resonators, it would excite that resonator. Those resonators which are 
excited correspond to factors of the target number. Resonators would be 
excited in conjugate pairs since each divisor has a corresponding quotient. 

Such a set of resonators could be constructed to suit any type of wave in any 
frequency range of choice. This is merely a spectrum analysis approach to 
factorisation in which the set of integers specifies the spectrum. An 
equivalent approach to constructing a physical spectrum analyser is to apply 
standard Fourier signal processing techniques. 

In practice, there are limitations on the precision, resolution and 
selectiveness of resonators whether implemented physically or as a software 
equivalent. Spectrum analysers are commonplace, and they have strong 
conceptual appeal. However it is worth bearing in mind that the factorisation 
spectrum is expected to be sparse. The most difficult factorisations would 
have only one pair of factors, and the spectrum could be analysed equally 
well by a progressive scan approach. 

Factorisation by a Frequency Analogue 

Instead of defining the wavelength of the fundamental in terms of the target 
number, it is possible to regard the target number as a frequency which is 
some harmonic of a divisor frequency. This analogue presents an attractive 
possibility illustrated in Figure 13. 

The divisor is generated as a square wave by a variable frequency oscillator. 
All the harmonics representing odd multiples of the divisor are then 
contained in the square wave signal. If one of these harmonics is equal to the 
target number frequency, then a factorisation has been identified. 

The higher the harmonic, the lower the amplitude and in order to establish 
whether a matching harmonic exists, the divisor signal is passed through a 
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very selective filter. This removes all harmonics except those which have 
frequencies in the vicinity of the target number frequency. 

The output from the filter is compared to the target frequency by means of a 
phase comparator. This generates an error signal which is used to control the 
frequency of the divisor oscillator. When one of the harmonics matches the 
target frequency exactly, the divisor frequency will lock. The frequency of 
the divisor oscillator can be measured precisely by a digital frequency 
meter. 



A phase locked loop factoriser. 

Figure 14. 

The frequency based analogue is in many ways simpler and more direct than 
the spectrum analysis of wavelength. 

Factorisation by Phase Equality 

When two waves coincide in space, they interfere with one another. When 
the waves have identical frequencies and a fixed phase relationship, 
throughout the region of space being observed, they give rise to standing 
wave patterns. Everyone is familiar with the phenomenon of interference of 
waves. The patterns created by waves as they are reflected back and forth 
from the banks of ponds are so commonplace as to be mundane. Yet they 
display interference as graphically as one could wish. Perhaps the most 
powerful and highly developed tool available to science is the 
interferometer. 

Figure 13 illustrated that factorisations can be represented by arcs of a circle 
lying on the perpendicular bisector of a line, the length of which represents 
2 Jt . This simply means that both ends of the line are equidistant from the 
centre. Thus, if a source of waves was placed at any point on the 
perpendicular bisector it would result in the waves arriving at both ends of 
the line in phase. 
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A factorisation requires that the radius of the circles have an integer value. 
This can be interpreted to mean that an exact number of waves can be 
counted between the source and the ends of the line. 

Factorisation also requires that the distance from the source to the midpoint 
of the line representing 2j~Tis also an integer value. So a different, but exact 
number of waves must match this distance too. 

An appropriate scale factor can be chosen so that the length 2fTis related 
to the wavelength of the source of waves. Under these conditions, placing 
the source at a factorisation point will result in the signal at the centre of the 
line, being in phase with that at either end of the line. 

Conversely, if a signal is emitted from each end of the line and from the 
centre of the line so that all three signals are emitted in phase, they will 
again be in phase at any factorisation point on the perpendicular bisector. 

By a symmetry consideration, it is possible to dispense with one of the 
signals from the line end. This leads to the arrangement shown in Figure 15. 



factorisation points. 

Figure 15. 
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Experiment 1. 


The principle of operation can be demonstrated readily using sound. A 
loudspeaker can be fed with a signal which is a pure tone and placed on the 
axis X-X. Any convenient sound frequency can be used. The speed of sound 
in air is approximately c=330»w _1 . Taking a signal of f=\.5KHz the 
wavelength will be: 



, 330 

1500 


0.22 


/ = 0.22 m 


If the number 105 is used as in Figure 12, the distance between SI and S2 
should be: 


s = ijm 

S = 0.22 x 10.247 = 2.254m 

A pair of microphones may be set to the above separation, at the locations 
SI and S2. The loudspeaker can be moved along the line X-X which is at 
90° to the microphones and in line with one of them. The strength of the 
signal at the microphones will vary and will decline overall as the 
loudspeaker is moved further from the microphones. 

However, the signals from the microphones will be found to be in phase 
with each other (and with the sound emitted at the loudspeaker), when the 
loudspeaker is at distances corresponding to the factorisations. These are: 
0.44m, 0.88m, and 1.76m respectively for the above conditions. 

Evidently, the distance between SI and S2 could be fixed as an alternative 
and the frequency of the sound altered to achieve the correct scaling. 

For best results, the experiment should be set up where there are not likely 
to be any strong reflections of the sound, perhaps in the open air. The size of 
the loudspeaker should be as small as possible so that it represents a point 
source. Ensure that both microphones are connected with the same polarity, 
and the sound used should be as pme a tone as possible. If the loudspeaker 
is highly directional, it may be preferable to place the loudspeaker on its 
back. This may help to produce a more even sound field in the plane of the 
microphones. Since the microphones are not the same distance from the 
loudspeaker, there will be a difference in signal strength between the 
microphones. 
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General purpose recording software can be used to examine the microphone 
signals. Phase is the important characteristic to observe. This is unaffected 
by any difference in loudness between the signals at the two microphones. 

A pure tone is helpful in understanding the operation of the factoriser and so 
preferable in the experiment described. However, since phase is the 
important characteristic being observed, a pulse train can be substituted for 
the pure tone. The pulse repetition frequency is then the surrogate for 
wavelength. The time interval between arrival of pulses becomes the means 
of phase comparison. 

Although carried out using sound, the experiment is very similar to a 
Young's Slits configuration which was devised by Thomas Young (1773 - 
1829), and gave the first effective demonstration of the wave nature of light. 

In Young's experiment, the plane of observation may be any plane parallel 
to the plane containing the slits. This means that the amplitudes of the 
signals from each of the slits are virtually identical. Thus, when the signals 
were in phase, the brightness would be double that from a single slit, and 
when the signals were out of phase, cancellation would be essentially 
complete. 

This an important distinction between Young's experiment and the 
configuration described here. For the purpose of factorisation the line of 
observation is coincident with one or other of the slits. The observation line 
is also normal to the plane containing the slits. This may be described most 
appropriately as an Axial Equiphase Factoriser. The distance of any 
factorisation point to each of the slits is different. Consequently signals from 
the slits at any factorisation point are not of identical intensity. Cancellation 
and addition of the signals does not produce complete extinction or twice 
the brightness. 

The diagram of Figure 15 is a two dimensional representation. The 
factorisation axis is nevertheless an axis of symmetry. A configuration with 
a central hole on the axis for SI and an annular aperture replacing the 
second slit S2 could be used. This latter arrangement is similar in 
construction to that used for a Fresnel lens. 

The type of wave or pulse train used in an Axial Equiphase Factoriser is a 
matter of choice, and a radio frequency based system could be employed. 
Very large radio astronomy arrays are in current use. There is in fact some 
similarity between the factorisation arrangement and some radio navigation 
systems. 
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The physical size of a factoriser based on interferometry is determined by 
the wavelength used, making an optical system perhaps the natural first 
choice in the design of a factoriser intended to handle large numbers. 
Tuneable lasers would then be the natural first choice as a light source. The 
path from SI to a factorisation point and that from to S2 to a factorisation 
point represent exact wavelength multiples and so satisfy the requirements 
for a laser cavity. 

There are many refinements which would be incorporated in a precision 
analogue factoriser based on the interference principle. The technical 
considerations are well known amongst opto-mechanical engineers. At the 
time of writing, the largest optical interferometers are those intended for the 
detection of gravitational waves. Although such waves may not exist, or may 
be undetectable by direct methods, the technology used in the 
interferometers demonstrates beyond doubt the degree of refinement which 
can be achieved with this type of instrument. 

Perhaps the greatest attraction of an Axial Equiphase Factoriser is that there 
is potential for the solution to become available in a time scale essentially 
unrelated to the complexity of the problem. 

It can be seen that the analogues described may be implemented in 
technologies ranging from the mechanical to the quantum mechanical, but 
even with an Axial Equiphase Factoriser it would not be reasonable to 
expect an analogue approach to determine a factorisation with absolute 
precision in any but the simplest cases. 

However, an analogue might serve to indicate those regions of the integer 
spectrum where a factorisation is most likely to be found. Numerical 
techniques would then be applied in those regions to identify the factors 
exactly. 
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Computer Architecture and Software Considerations 

The investigation described was carried out on a standard microcomputer. 
The largest integer handled by the compiler employed was a signed 64 bit 
type. The limit on size of integer makes practical comparison between the 
different algorithms used somewhat unsatisfactory. 

The processor registers on the computer used do not handle 64 Bit numbers 
explicitly. The memory structure also requires additional address cycles to 
transfer the extra bytes of data associated with numbers of this size. 

There is no doubt that manufacturers have devoted considerable attention to 
the optimisation of multiplication and division both from a hardware and 
software perspective. It is to be expected that for these operations at least, 
the compiler would make use of the full set of processor registers and cache 
memory when manipulating 64 Bit numbers rather than rely on memory 
cycles to transfer the data as required. Indeed, multiplying circuitry which 
can carry out a floating point multiply operation in a single cycle has been 
available for many years. Ever increasing electronic integration has 
extended the capacity of such devices. At the time of writing, commercially 
available Digital Signal Processors (DSPs) carry out a 32 bit multiply in a 
single cycle. 

Addition, subtraction, and comparison of numbers have traditionally been 
regarded as fast operations. It seems unlikely that these operations will have 
benefited from the same commitment to optimisation. Consequently it is 
likely that the comparisons made herein are flawed in favour of direct 
division approaches. 

It should also be remembered that the unavoidably large number of calls to 
library subroutines brings with it the overhead of context switching. 
Embedded code is preferable. It has also been noted that library subroutines, 
by their very nature, cannot be optimised to the needs of algorithms such as 
those described in this volume. 

Without becoming embroiled in direct coding, 64 bit numbers represented 
the limit of the comparisons available. The machine architecture of home 
computers was never intended for handling integers of arbitrary size. In 
order to take algorithm development further the natural progression would 
be construction of a dedicated factoriser card controlled by a microcomputer 
or DSP. 
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Number Bases 


All of the examples which have been considered are expressed to base 10. 
General purpose digital computers use binary, which confers tremendous 
flexibility in use. The algorithms described have naturally been programmed 
in the binary context. The problem of factorisation can be addressed in any 
number base, but that does not mean that the performance will be the same 
for all number bases. Furthermore, the relative problem complexity in 
different number bases is not independent of the problem size. It is entirely 
feasible to construct a factoriser which operates in base 10. 

Considerations for a Purpose Built Factoriser. 

Addition and subtraction are generally faster operations than either of 
division or multiplication. The speed of large scale integration (LSI) logic 
circuits commercially available at present is typically up to about 125 MHz, 
though some circuits can operate much faster than this. The clock rate above 
refers to readily available logic circuits which can be incorporated into the 
core of a special purpose factoring engine. Within such a computing engine, 
the least significant data bits must be able to change much more rapidly than 
their higher order counterparts. Consequently, the overall speed of any 
design can be increased through the use of carefully designed special 
purpose prescaler circuits to handle the most rapidly changing data bits. 

The machine architecture required for the decrementing algorithm is 
straightforward, and is ideally suited for parallel computation. Similarly, if it 
was desired to construct a moderate sized sieve of Eratosthenes the use of 
wired OR modulo N counters is a practical proposition. 

In the case of combined algorithms, the flexibility of central processor 
supervision of dedicated circuitry has an undeniable attraction. The 
algorithms themselves are relatively compact whereas the numbers being 
manipulated can be arbitrarily large. This suggests that some form of 
Harvard architecture might be appropriate. The separate memory and bus 
structure for program control offers the designer flexibility. A width of data 
bus appropriate to the problem size is more readily accommodated in a 
Harvard architecture than in a monolithic memory structure. It is 
conceivable that an ordinary binary mode processor could supervise 
calculations carried out directly in other number bases without concern for 
conversion between bases. 
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Epilogue 

The factorising problem has its roots in the infinity of possible factors for 
any real number. Requiring integer solutions is an arbitrary restriction for 
which there does not appear to be any convenient analytical technique. In 
the absence of a purpose built factorising engine, recourse is made to the 
microcomputer, and while this imposes may limitations, it provides many 
interesting opportunities for investigation. 

In practice, for numbers up to about 10 12 , the performance of direct division 
surpassed that of the decrementing algorithm described here. However, for 
numbers in the region of 10 18 , the situation was decidedly reversed with the 
decrementing algorithm described running between 5 and 17 times faster 
than a simple direct division algorithm, depending on the actual composition 
of the number being factorised. 

As expected, the algorithms described are most suitable for factorising large 
numbers, with the advantages becoming more pronounced as the size of the 
integer is increased. The comparison cited above was not absolute. In the 
case of direct division, it was necessary to employ 64 bit variables. The 
decrementing algorithm was able to take advantage of smaller variable types 
yet retain the capacity to factorise large numbers. The computational 
overhead resulting from the use of 64 bit variables has already been noted. 

Nevertheless, the variation in performance as the size of integer was 
increased was quite clear. The distinct shift in advantage from direct 
division to the decrementing algorithm was independent of the variable 
types used and would have reflected the characteristics of the algorithms 
themselves. 

Comparison could not take account of the degree to which hardware and 
compilers have been optimised for division and multiplication, and the 
software implementation of the algorithms discussed here is thought to be at 
a decided disadvantage to direct division. 

No matter how powerful computers may become, it will always be possible 
to specify a factorisation problem which exceeds the available capacity of 
single cycle multipliers. For an arbitrary length integer, recourse must 
ultimately be made to traditional software solutions which handle 
multiplication and division in a piecewise manner over many cycles. 

In contrast, the hardware for a decrementing algorithm approach can be 
extended cheaply and, for all practical purposes, indefinitely. The 
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computational burden is proportional to JT, which although feeble in 
comparison to a direct analytical solution is quite respectable. 

The approaches based on quadratic formulations have potential application 
in respect of numbers much larger than those examined in this study. Their 
incorporation, together with finite difference enhancements would serve to 
reduce the proportionality coefficient of Jr referred to above. 

Overall, the performance of the algorithms was gratifying and suggests 
considerable potential for any purpose built factorising engine which 
embodies techniques similar to those described in this volume. 

The prospect of an almost real time hybrid factoriser based on 
interferometry is really quite interesting. The history of Science and 
Mathematics persistently reveals that Nature has convenient means of 
answering our questions even before mankind has formulated the problems. 

Judicious application of rational thought and careful investigation confirms 
that factors are ultimately condemned to reveal their true values. 
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Appendices 

The appendices contain the core code of the programs which are written in 
Pascal. The user interface declarations have been omitted since these are a 
matter of style and primarily reflect the particular development environment 
employed. 

The programs were written in order to confirm the validity of the algorithms 
and were deliberately given a monolithic structure. No particular attention 
has been paid to exception handling. Certain parts might have been 
structured as procedures without impairing performance, but maintainability 
was never a concern and readability will always be a matter of personal 
taste. 

A compiled program running on a general purpose commercial operating 
system could never be considered a serious contender by anyone who 
wished to stretch the limits of factorisation methods. 

The programs will hopefully serve as an effective introduction to the 
techniques described. 
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Appendix 1 - Basic Decrementing Algorithm 


This is the original decrementing algorithm written on 26th March 2007. It 
has no optimisation. The program serves as a reference point to confirm the 
correctness of the faster algorithms. This program is designed to find a pair 
of factors for any given number. The factors may be either composite 
numbers or primes. Repeated application of the algorithm can be used to 
find the prime decomposition if desired. 

VAR 

Targetnumber, Divisor, Quotient, Remainder: Longint; 

(* The number to be factorised is entered. The largest number which 
Longint can accommodate is about 10 digits. The type "COMP" is a larger 
integer with about 20 digits, but is classed as a real and is not compatible 
with certain integer operations.*) 

Begin 

Targetnumber := Strtoint(Editl.Text); 

Divisor := Trunc(Sqrt(Targetnumber)); 

Quotient := Targetnumber DIV Divisor; 

Remainder := Targetnumber MOD Divisor; 

While (Remainder <> 0) Do 
Begin 

Divisor := Divisor - 1; 

Remainder := Remainder + Quotient; 

While (Remainder >= Divisor) Do 
Begin 

(* The number of times that this loop is traversed increases as the 
divisor becomes smaller and the quotient grows in size. Initially 
the loop is completed once. If the target number is a prime, then 
the loop is traversed about target/2 times. The program quickly 
returns 9 (3x3) and 220850491 (prime) as factors of the number 
1987654419. However, the number 1987654421 is prime, and the 
program takes noticeably longer. *) 

Remainder := Remainder - Divisor; 

Quotient := Quotient + 1; 

End; (* While remainder > Divisor *) 

End; (* While remainder not zero *) 

End; 
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Appendix 2 - Improved Decrementing Algorithm 


This is the improved decrementing algorithm using the multiple and 
remainder. The algorithm handles numbers up to about 10 A 19. All the 
variables within the main search loop are defined as Longint. This is 32 bits 
wide. Defining these as COMP increases the capacity of the program in 
principle. However the computational burden incurred reduces program 
speed. 

Overflow of the longint variables in the main loop can cause the Remainder 
to go negative causing program failure in the case of primes > 10 A 12. 

VAR 

Targettext: String; 

Smallprime : Array[1..7] of Integer; 

Quotient, Targetnumber, Primeprod, Tempnum, Primemod : Comp; 

Divisor, Remainder, Diffrence, Primetest: Longint; 

Multiple, Twomult: Longint; 

Lastone, Digittest, I, Cnt: Integer; 

Factored : Boolean; 

(* NOTE For a hard wired implementation of this algorithm, word length is 
essentially irrelevant and would not compromise computational speed. 
Furthermore, some of the loop overheads would be eliminated and several 
operations carried out in parallel. Even without going to the trouble of hard 
wiring, direct coding could make full use of the registers and eliminate many 
memory accesses and subroutine calls inevitable with compiled code. *) 

Begin 

Factored := FALSE; 

Quotient := 0; 

Divisor := 1; 

Lapsed := 0; 

Targettext := Editl.text; 

Targetnumber ;= Strtofloat(Targettext); 

(* The following tests for divisibility can be carried out by direct division.*) 

Lastone ;= Length(Targettext); 

Digittest := Strtoint(Copy(Targettext,Lastone,l)); 
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If (Digittest = 5) Then 
begin 

(* Divisible by 5*) 

Factored := True; 

Divisor := 5; 
end 
Else 
begin 

If (NOT odd(Digittest)) Then 
begin 

(* Divisible by 2 *) 

Factored := True; 

Divisor := 2; 
end 
Else 
Begin 

For Cnt := 1 to (Lastone - 1) Do 
Begin 

Digittest := Digittest + Strtoint(Copy(Targettext,Cnt,l)); 

End; (* Adding the digits *) 

If ((Digittest MOD 3) = 0 ) Then 
Begin 

(* Divisible by 3 *) 

Factored := True; 

Divisor := 3; 

End 

Else 

Begin 

(*-Divide target number by product of low primes-*) 

(* The product of primes 7,11,13,17,19,23 = 7436429. If non zero, the 
remainder from division by this number into the target number is then 
divided by each of the primes in turn to ensure that none of them divide 
the target number. With a conventional computer, division time depends 
on the size of the numbers. It should be more efficient to reduce the size of 
the division by the above approach than to divide each prime in turn into 
the target number. In practice it makes no significant difference until a 
very large number of such direct divisions are being replaced. *) 

Smallprime[l] := 7; 

Smallprime[2] := 11; 
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Smallprime[3] := 13; 

Smallprime[4] := 17; 

Smallprime[5] := 19; 

Smallprime[6] := 23; 

Smallprime[7] := 29; 

Primeprod := 215656441; 

(* Primeprod is the product of the small primes. It is defined as COMP 
to ensure compatibility with targetnumber. *) 

Primemod := Int(Targetnumber / Primeprod); 

Primetest := Trunc(Int(Targetnumber - (Primemod * Primeprod))); 

If (Primetest = 0) Then 
Begin 

Factored := TRUE; 

Divisor := 7; 

(* Any mix of the small primes could be used instead. 7 has been 
chosen as the smallest divisor. *) 

End 
Else 
Begin 
i := 1; 

While (Not(Factored) AND (i <= 7)) Do 
Begin 

Divisor := Smallprime[i]; 

If (Primetest MOD Divisor = 0) Then 
Begin 

Factored := TRUE; 

(* The divisor has a value assigned already. *) 

End 
Else 
Begin 
i :=i+ 1; 

End; (* If primetest *) 

End; (* While not factored *) 

End; (* If primetest = 0 *) 

End; (* For cnt := 1 *) 

End; (* If not odd digittest. *) 

End; (* If digittest = 5 *) 
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(*-This completes the trivial and low prime tests.-*) 

If (Factored) Then 
Begin 

Quotient := Int(Targetnumber / Divisor); 

(* Quotient is left as type COMP since it is not used in the main 
program loops. *) 

Remainder := 0; 

(* Note that this may be only one of many factorisations possible. *) 

End 

Else 

Begin 

(* no easy factorisation has been found so initialise main search. *) 

Divisor := Trunc(Sqrt(Targetnumber)); 

(* The square root is taken as the starting point for simplicity. The 
search can be started from any divisor, making the algorithm suitable 
for parallel implementations. *) 

Diffrence := 0; 

Tempnum := Divisor; 

(* Tempnum (COMP) is used to avoid overflow when the divisor is 
squared. *) 

Remainder := Trunc(Targetnumber -(Tempnum * Tempnum)); 

While (Remainder >= Divisor) Do 
Begin 

(* This loop takes care of the approximations incurred by rounding 
down the square root and notionally equating the Divisor and Quotient. *) 
Remainder := Remainder - Divisor; 

Diffrence := Diffrence + 1; 

End; (* While remainder *) 

If Divisor <= 30 Then 
begin 

(* All divisors up to 29 have been tested already. The next valid 
test is for divisor = 31. i.e targetnumber >=961. Consequently, a 
divisor <= 30 implies that the number is prime. *) 

Divisor := 1; 

End; (* If divisor <= 30 *) 
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(*-Main Search Loop-*) 

Multiple := 1; 

Twomult := 2; 

(* Multiple is quotient MOD divisor, and is initially unity. As the divisor 
diminishes. Multiple increases. Multiple is actually the number of times 
the remainder / divisor loop would be traversed in the simple algorithm. 
The smaller the divisor, the greater the number of loops saved.*) 

While ((Remainder > 0) AND (Divisor > 30)) Do 
Begin 

Divisor := Divisor - 1; 

Remainder := Remainder + Diffrence + Multiple; 

(* The remainder is affected by the old difference *) 

Diffrence := Diffrence + Twomult; 

(* The new difference is used next time round. *) 

While ((Diffrence >= Divisor) AND (Remainder >= Divisor)) Do 
Begin 

(* The difference between the size of the quotient and the divisor 
has exceeded the size of the divisor. So Quotient DIV Divisor has 
increased by at least unity. *) 

Diffrence := Diffrence - Divisor +1; 

(* The addition of unity corresponds to the increment in the quotient 
required every time an extra subtraction of the divisor is required. *) 
Multiple := Multiple + 1; 

Twomult := Multiple * 2; 

Remainder := Remainder - Divisor ; 

End; (* While Diffrence *) 

While (Remainder >= Divisor) Do 
Begin 

(* The number of times that this loop is traversed changes as the 
divisor becomes smaller. *) 

Remainder := Remainder - Divisor; 

Diffrence := Diffrence + 1; 

End; (* While Remainder >= Divisor *) 

End; (* While Remainder > 0 and divisor > 30 *) 

(*-End of Search Loop-*) 

(*-Finish off-*) 
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If (Remainder = 0) Then 
begin 

(* Normal termination of algorithm.*) 

Quotient := Int(Targetnumber / Divisor); 

Tempnum ;= Divisor; 

Remainder := Trunc(Targetnumber -(Tempnum * Quotient)); 

(* Calculating remainder confirms that the division algorithm was 
correct.*) 

End; (* If remainder = 0 *) 

If (Divisor <= 30) Then 
Begin 

(* All prime numbers between 1 and 28 have been tested. Thus the 
targetnumber 
is a prime. *) 

Divisor := 1; 

Quotient := Int(Targetnumber); 

Remainder := 0; 

End; 

End; (* Else part of If factored *) 

If (Remainder < 0) Then 
Begin 

MessageDlg('Factorisation process failed'.MtWanting,[Mbok],0); 
(* Factorisation fails through overflow causing the remainder to 
become negative. *) 

End; 

End; 

end. 
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Appendix 3 - First Quadratic Algorithm 


VAR 

Targettext: String; 

Factored : Boolean; 

Rootl, Root2, Targetnumber, Mu, Twomu, Delta, Lambda: Comp; 
Deltasqr, Lambdasqr, TwoLambdaplus, Upperlim : Comp; 

FourmT, FourT, Roottgt: Comp; 

Deltal :Real; 

Begin 

Factored := FALSE; 

Rootl := 1; 

Mu := 1; 

Targettext := Editl.text; 

Targetnumber := Strtofloat(Targettext); 

Roottgt := Sqrt(Targetnumber); 

(*-Trivial Divisibility Tests Inserted Here See Appendix 2-*) 

If (Factored) Then 
Begin 

Root2 := Int(Targetnumber / Rootl); 

(* Note that this may be only one of many factorisations possible. *) 
End 
Else 
Begin 

(* no easy factorisation has been found so initialise main search. *) 
Lambda := Int(Sqrt(Targetnumber)); 

(* The square root is taken as the starting point for simplicity. *) 

If ((Targetnumber - (Lambda * Lambda)) = 0 ) Then 
Begin 

(* The targetnumber is a perfect square, so there is nothing further 
to be done*) 

Rootl := Lambda; 

Root2 := Lambda; 

End 

Else 

Begin 

(* This is the full search loop. *) 
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Upperlim := Targetnumber/(Mu + 1); 

(* Delta is the difference between Mu*Divisor and the quotient. When 
Delta exceeds the value of the Divisor, Mu should be increased by one 
and Delta reset. The point at which Delta exceeds the divisor decreases 
as the divisor decreases. In general, this occurs at Sqrt(T/(Mu+l)). Since 
Deltasqr is readily available throughout, the upperlimit on it is used as 
a breakpoint instead. Thus Upperlim = T/(Mu +1). *) 

FourT := 4 * Targetnumber; 

FourmT := FourT; 

(* Mu is initially unity. *) 

Lambda := Int(Sqrt(FourmT)) +1; 

Lambdasqr := Lambda * Lambda; 

Deltasqr := Lambdasqr - FourmT; 

TwoLambdaplus := (2 * Lambda) + 1 ; 

Deltal ;= Sqrt(Deltasqr); 

Delta := Int(Deltal); 

While ((Deltal - Delta) > 0 ) Do 

(* The search ends when Deltasqr is a perfect square. *) 

Begin 

If (Deltasqr > Upperlim) Then 
Begin 

(* Increase mu by unity and adjust other variables.*) 

Mu ;= Mu + 1; 

FourmT := FourmT + FourT; 

Lambda := Int(Sqrt(FourmT)) +1; 

(* Calculation of Lambda is probably the most computationally 
intensive operation in this initialisation section. *) 
TwoLambdaplus := 2*Lambda + 1; 

(* Multiplication by 2 should be a simple shift operation.*) 
Lambdasqr := Lambda * Lambda; 

Deltasqr := Lambdasqr - FourmT; 

Upperlim := Int(Targetnumber/(Mu + 1)) + 1; 

End 

Else 

Begin 

Deltasqr := Deltasqr + TwoLambdaplus; 

TwoLambdaplus := TwoLambdaplus + 2; 

End; (* If delta <= upperlim *) 

Deltal := Sqrt(Deltasqr); 


73 



Delta := Int(Deltal); 

End; (* While Deltasqr not a perfect square. *) 
Lambda := Sqrt(Deltasqr + FourmT); 

Twomu := 2*Mu; 

Rootl := (-Delta + Lambdaj/Twomu; 

Root2 := Targetnumber / Rootl; 

End; (* Else part of If Targe tnumber etc. *) 

End; (* Else part of If Factored. *) 

End; 

end. 
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Appendix 4 - Second Quadratic Algorithm 

Quadratic factoring. This version implements repeated application of the 
quadratic formula to locate possible convergence points. 

VAR 

Targettext: String; 

Factored : Boolean; 

Rootl, Root2, Targetnumber, Mu, Twomu, Fourmu: Comp; 

Lambda, Lowerlim Divisor, Quotient, Remainder, Deltaplusl: COmp; 
Deltaplussqr, Possroot, FourmumrgnO, Diffrence, Rootdiffl: Comp; 
Rootdiff2, Lastdiff, Lastroot, DivDeltaRootDif: Comp; 

Discrimrt :Real; 

Begin 

Factored := FALSE; 

Rootl := 1; 

Mu := 1; 

Targettext := Edit 1.text; 

Targetnumber := Strtofloat(Targettext); 

(*-Trivial Tests for Divisibility are inserted here-*) 

If (Factored) Then 
Begin 

Root2 := Int(Targetnumber / Rootl); 

(* Note that this may be only one of many factorisations possible. *) 
End 
Else 
Begin 

(* no easy factorisation has been found so initialise main search. *) 
Lambda := Int(Sqrt(Targetnumber)); 

(* The square root is taken as the starting point for simplicity. The 
search can be started from any divisor, making the algorithm suitable 
for parallel implementations. *) 

If ((Targetnumber - (Lambda * Lambda)) = 0 ) Then 
Begin 

(* The targetnumber is a perfect square, so there is nothing further 
to be done*) 

Rootl := Lambda; 

Root2 := Lambda; 
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End 

Else 

Begin 

Divisor := Lambda; 

Lowerlim := Int(Lambda/Sqrt(Mu +1)); 

Quotient := Int(Targetnumber/Divisor); 

Remainder := Targetnumber -(Divisor * Quotient); 

(* Various initialisations... *) 

Twomu := 2 * Mu; 

Fourmu := 2 * Twomu; 

Possroot := 51; 

Lastroot := 0; 

Lastdiff := 0; 

While ((Remainder <> 0) AND (Divisor > 1) AND (Possroot >50)) Do 
Begin 

(* Mu is initially unity. *) 

(* Calculate a new possible root.... The quadratic formula as applied to 
this problem. This is a computationally onerous part of the program. *) 

Deltaplusl ;= Quotient - (Mu * Divisor) + 1; 

Deltaplussqr := Deltaplusl * Deltaplusl; 

FourmumrgnO ;= FourMu * (Divisor - Remainder); 

Discrimrt := Sqrt(Deltaplussqr + FourmuMrgnO); 

Possroot := Int((-Deltaplusl + Discrimrt)/Twomu); 

(* The separation of the possible roots is a fairly well behaved function 
it should be practical to replace application of the quadratic formula by 
a finite difference extrapolation without endangering the factorisation 
process. *) 

Rootdiffl := Possroot - Lastroot; 

Lastroot := Possroot; 

Rootdiff2 := Rootdiffl - Lastdiff; 

Lastdiff := Rootdiffl; 

(* Update the divisor etc.*) 
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Divisor := Divisor - Possroot; 

Quotient := Quotient + Possroot; 

Remainder := Targetnumber -(Divisor * Quotient); 

(* After the first few steps, the gap between successive possible roots 
declines. When it becomes less than unity the algorithm fails. *) 

While ((Remainder > Divisor) AND (Divisor >1)) Do 
Begin 

(* This is simply to ensure that the value of the divisor and quotient 
produce a remainder in the correct range. Recall that the divisor is the 
control variable and the value of n determines the amount by which the 
divisor is reduced. It is a matter of convenience to increase the quotient 
by n. This may not be the correct amount but will always be less than or 
at best equal to the correct quotient. Consequently the remainder may 
turn out to be greater than the divisor. *) 

Remainder := Remainder - Divisor; 

Quotient := Quotient + 1; 

End; (* While remainder > Divisor.... *) 

While ((Remainder < Divisor) AND (Divisor > 1)) Do 
Begin 

(* The value of n has brought us close to where the remainder will 
exceed or equal the divisor. It should be within unity of the transition, but 
may not be on account of the rounding processes which have been 
applied. This decrementing loop adjusts the divisor, quotient and 
remainder until the remainder exceeds or equals the divisor. *) 

Divisor := Divisor - 1; 

Quotient := Quotient + 1; 

Diffrence := Quotient - Divisor - 1; 

Remainder := Remainder + Diffrence; 

End; (* While remainder < Divisor.... *) 

(* Once the remainder has exceeded the divisor, the quotient has to be 
incremented while leaving the divisor the same. The remainder is also 
adjusted.*) 

Quotient := Quotient + 1; 

Remainder := Remainder - Divisor; 

If (Divisor < Lowerlim) Then 
Begin 

Mu := Mu + 1; 
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Twomu := 2 * Mu; 

Fourmu := 2 * Twomu; 

Lowerlim := Int(Lambda/Sqrt(Mu +1)); 

End; (* If divisor < Lowerlim *) 

Loopcnt := Loopcnt +1; 

DivDeltaRootDif := Int(Divisor/(Deltaplusl - l))-Possroot; 

End; (* While possroot > 50 *) 

Rootl := Divisor; 

Root2 := Quotient; 

End; (* Else part of If Targetnumber etc. *) 

End; (* Else part of If Factored. *) 

End; 

end. 
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